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FOREWORD 


This report summarizes all work performed by Southwest 
Research Institute during the past two years, under Contract No. 
NAS8-25919, "Analysis of Propellant Feedline Dynamics. 11 This study 
was performed for the George C. Marshall Space Flight Center of 
the National Aeronautics and Space Administration, and was adminis- 
tered technically by Messrs. Larry Kiefling and Gary Muller of the 
Aero-Astrodynamics Laboratory. 
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I. INTRODUCTION 


I. 1. Background, and Objectives 

At least three classes of feed system instabilities have been 
identified for liquid propellant rockets 3, 4->' < in eac h case, the 
dynamic behavior of the feed line plays a major role in the instability. 
Probably the best known problem^-* ^ involves a closed-loop coupling 
of the vehicle longitudinal structural modes, the feed system and the 
engine. A second problem^ showed up on Saturn flight AS- 5 04 when 
the S-Il center engine displayed severe pressure and thrust oscilla- 
tions rather late in the stage burn. In this case, the problem was 
localized, involving a closed-loop instability of the engine support 
structure, the feed system, and the engine. A third type of problem 
which can exist^ involves only the feed system, the inducer -pump 
combination and the engine. This instability is traceable directly to 
the inducer dynamics, with the feed system and engine providing an 
impedance loading which governs the resultant oscillation frequency. 

All of the above instabilities pose a threat to the success of a 
mission. For simulation purposes, it is readily recognized that ade- 
quate modeling procedures are needed for each portion of the vehicle 
involved in a given instability. The objective of this study was to de- 
velop an analytical method, and its corresponding computer program, 
to allow the study of disturbances of liquid propellants in typical engine 
feedline systems. This method was to include (1) the effect of steady 
flow, (2) the influence of distributed compliances and axial wall stiff- 
ness, (3) the effects of local compliances, such as vapor bubbles, and 
(4) various factors causing coupled responses between the line and 
structure, i. e. , bends, mounting stiffness, forced changes in line 
length, and bellows and PVC couplings. 

An additional objective of this program was to construct a 
basic test apparatus to simulate a feedline system, consisting of an 
upper reservoir, a feedline, a terminational resistance (impedance), 
a lower reservoir, and a circulation pump; this- test facility was to 
be used to experimentally verify many of the various effects dis- 
cussed above on the frequency response of a feedline. In general, 
the experimental program was to provide a check-out or verifica- 
tion of the analytical model and computer program. 


Superscript numbers in the text refer to References presented on 
page 148 of this report. 
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In all respects, this objective has been satisfied. This report 
summarizes the modeling results for each item defined above, and 
presents a generalized computer program which allows the user to 
easily deter min e dynamic .performance of a—line containing any number 
of elemental components. Use of this computer program is illustrated 
by way of a detailed set of instructions contained in Appendix D, plus 
a presentation of several example problems. 


I. 2. General Feedline Problem 


The problem of modeling a given feedline can best be illus- 
trated by an example. Figure 1 shows a hypothetical case involving 
a propellant tank, a feedline, a pump and an engine. In general, 
each element of the model is elastically supported and can experi- 
ence some motion relative to a frame of reference on the vehicle. 
These accelerations, plus pressure oscillations in the engine com- 
bustion chamber and cavitating inducer instabilities, represent dis- 
turbance inputs to the feed system which can cause flow and pres- 
sure oscillations. These flow disturbances may, in turn, couple 
back with the engine and/or structure to produce a large amplitude, 
closed- loop instability. The modeling problem, with respect to the 
feed system, is to describe analytically the dynamic pressure and 
flow at an arbitrary point in the system, with proper treatment of 
the end conditions, the acceleration inputs and the flexible supports. 

Generally, a feedline may be regarded as a series of elements, 
such as lengths of line, bend's, bellows PVC coupling, etc. There- 
fore, it is necessary to describe each element individually, and then 
define means for mathematically combining these elements so as to 
produce the proper overall system mode. In this report, each ele- 
ment is basically defined by a four-terminal representation in the 
Laplace domain. Overall system model formulation is also achieved 
in the Laplace or operator domain. Where problem solution is re- 
quired only in the frequency domain, this Laplace formulation is the 
one to use, and the present computer program is set up on this basis. 
For problem solution in the time domain, the Laplace operator for- 
mulation must be converted to an equivalent time domain formulation 
(rational approximate model). The procedure for accomplishing 
this conversion is described herein. 
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II. GENERAL LINE MODEL WITH MEAN FLOW 


II. 1 . Physical Picture and Requirements 

A 21 illustration of the physical picture of the flow behavior in 
a feedline is given in Figure 2. For the present, only the fluid dy- 
namic effects in a line, assumed cylindrical, straight, and with in- 
finitely stiff walls, will be considered. 


The line shown in Figure 2 is assumed to contain a turbulent- 
ly flowing viscous liquid. The mean vector velocity is represented 
by v 0 , which contains only the axial component, v 0 (r). Superim- 
posed on this mean flow are two dynamic velocity perturbation 
quantities, v and v t . The velocity component v represents the 
dynamic velocity perturbation of interest, and vt represents the 
turbulent velocity fluctuations. Similarly, the total pressure is the 
sum of a mean quantity, p 0 , a dynamic -perturbation part, p, and 
a turbulent part, p t . 


By substitution of the assumed summation forms for velocity 
into the Navier-Stokes equations^, and short-time averaging of the 
turbulence quantities, the following form results: 

l? +V “^ = V(VV)-VX(VXV)} + TDT (1) 

In addition, we may complete our description of the dynamic 
fluid behavior with a continuity equation: 


3p 

at 


+ Po 


V * V + v 0 


la 

Bx 


0 


( 2 ) 


and a liquid equation of state: 

dp = K (3) 

Equation (1) is a first-order form of the Navier-Stokes equations, 
including turbulent dissipation terms (TDT). Similarly, Equation 
(2) is a first-order continuity expression relating the dynamic velo- 
city perturbation to density perturbations (p 0 is the mean density). 

Examination of Equations {1) and (2) shows that the presence 
of a mean turbulent flow has two possible effects. First, the 



TURBULENT PROFILE 



v T ( r,x,t ) = vector total velocity 

V V v+ 7 t ' 

v 0 =time invariant mean flow 
v= dynamic perturbation of interest 
v t = turbulent fluctuations 

P T = Po + P + Pt 

2747 


Figure 2. General Physical Picture Of Flow Behavior In Feed Line 
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turbulent dissipation terms introduce damping in addition to the con- 
ventional "laminar" damping accounted for by the terms in brackets 
on the right side of Equation (1). Second, the mean flow introduces 
a "convective" effect represented by the terms v 0 dv/dx and 
Vo dp/dx in. Equations .(1) and (2). respectively... 


The requirements for a dynamic feedline model with a mean 
flow are that (1) the attenuation effects of the turbulent flow be pro- 
perly treated and (2) the possible convective influences be examined 
and included if necessary. The following sections treat these two 
topics. In general, the feedline model presented will be an exten- 
sion of an existing distributed parameter model which exactly ac- 
counts for laminar viscous effects, but negligible mean flow. This 
existing model, which is discussed in References 5 and 6, may be 
given in the following form: 

P 2 (s) = Pi (s) cosh yL - Z 0 A 1 Qj. (s) sinh yL 

Q 3 (s) = Qi(s) cosh yLi - Zo^APjfs) sinh yL (4) 


The laminar propagation operator takes the familiar form 


Ylam “ 



2Jj (?r Q ) 1~ a 

?r 0 J 0 (§r 0 ) J 


(5) 


where c is the phase velocity in the fluid, J 0 and Ji are zero and 
first-order Bessel functions of the first kind and 

= -s/v (6) 

The characteristic impedance is related to the propagation operator 
by 


Z c = p 0 c 0 s y/s (7) 

The magnitude of the additional attenuation due to turbulence 
is a function of the Reynold's number and the frequency range of 
interest; this aspect is discussed in Section IL 3. The relationship 
between mean flow velocity and the adiabatic speed of sound dictates 
the phase velocity to be used in Equation (5). Section II. 2 discusses 
the effect of mean flow convection. * 
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II. 2 Possible Convective Influences 


As indicated previously, the terms v 0 dv/Sx and v 0 dp/dx 
in Equations (1) and (2), respectively, represent the convective 
effect of the mean flow v 0 . The question to be answered at this 
point is: Are these terms of enough significance to be retained in 
the final model? Ideally we would' like to retain the convective terms 
at all times, but, from a practical standpoint, they increase con- 
siderably the complexity of the final solution. 


To judge the importance of the convective terms, it is con- 
venient to consider a nondissipative, one -dimensional form of Equa- 
tions (1) and (2) or 


dvx dvx 

at + Vo dx 


_L 

Po 


and 


ap 

at 


+ Po 



+ 



o 


( 8 ) 


( 9 ) 


Transformation of Equations (3), (8), and (9) into the Laplace domain, 
and the subsequent elimination of p and P, yields one equation in 
V x (V x is the transform of v x ), or 


<c 0 3 -v 0 2 ) 


d 3 V x 
dx s 


- 2s v 0 


s 3 V x 


= 0 


where 


( 10 ) 


c 0 = (ft/p 0 ) s 


Laplace variable 
isentropic speed of sound 


Vo 


= axial mean flow velocity 


Assuming a solution for Equation (10) of the form 

V x = Ce Yx (11) 

we find the characteristic equation to be 

( c 2 - v 3 ) Y 3 - 2s VoY - s 3 = 0 (12) 


whose roots or values for Y are 
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( 13 ) 


where c* 2 = cf - vf . The plus sign applies to waves traveling in 
the retrograde direction, and the negative sign to waves traveling in 
the mean flow direction. 


From these results, it is found that for propagation in the +x 
direction, with the mean flow, the phase velocity is greater than c 0 . 
To demonstrate this result, recall that the propagation operator is 
functionally equivalent to the Laplace variable divided by a velocity 
quantity. Substituting the equivalent expression for c* into Equa- 
tion (13), the following result is obtained: 


s (v 0 ± Co) 

(v 0 + Co) (Vo - Co) 


(14) 


For propagation in the mean flow direction. 



s 

Co + v 0 


(15) 


Therefore, the phase velocity for downstream propagation is greater 
than the adiabatic speed of sound. Similar reasoning verifies that 
the phase velocity for retrograde propagation is less than the adia- 
batic speed of sound. 

It is important to point out that for normal conditions of 
liquid propellants flowing in a feedline at velocities up to, say, 50 
fps, the phase velocity with mean flow is so minutely changed that 
this convective effect may be neglected. The reason for this is that 
the mean velocity is much smaller than the liquid sonic velocity, c 0 . 
The only condition under which this might change would be where 
sufficient entrained gas is present in the fluid to cause the new ef- 
fective c 0 value to be, say, an order of magnitude lower than for 
the pure liquid. Then the convective effect should be included. For 
example (see Section III. 1), if the level of entrained gas were to 
approach 0. 1 percent by weight of the total liquid mass, then the 
equivalent speed of sound of the mixture would be approximately 
100 fps. In this case, the speed of sound is only twice as large as 
the maximum mean flow velocity, and the phase velocity would be 
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100 ±50 fps, depending on the propagation direction. This example, 
however, is unrealistic because the gas concentration rarely ap- 
proaches a level where the mean flow velocity is but a small frac- 
tion of the mixture speed of sound. Therefore, for all practical 
purposes, the phase velocity can be taken to be the adiabatic speed 
of sound. 

II. 3 Effect of Turbulence on Damping and Phase Velocity 


It is known that the four terminal Laplace domain represen- 
- tation for pressure and flow in a line is extremely accurate when the 
flow is laminar and the propagation operator is of the following form: 


L = 


(*k\ 

' 2Ji(§r 0 ) “]' 

UJ 

§roJo(?ro) -1 


■* 


Ylam 


( 16 ) 


In reality, a turbulent flow condition exists for many propellant feed- 
line operating points, as indicated by Reynold's numbers of the order 
of 10 5 to 5 X 10 7 . Therefore, the effect of turbulence must be re- 
flected in the propagation operator. As shown in Figures 3 and 4 
(from’ Reference 7), the predominant effect of turbulence is to in- 
crease the spatial attenuation at low frequencies. At high frequen- 
cies, the laminar and turbulent attenuations coincide. There is some 
tendency for turbulence to reduce the phase velocity at very low fre- 
quencies and high Reynold's numbers. However, after a thorough 
investigation into existing feedline geometries, propellant properties, 
mean flow velocities and frequency range for structural-hydraulic 
coupling, it was concluded that the effect of turbulence on the phase 
velocity could be neglected. 


For the exact distributed parameter line model, the effect of 
turbulent attenuation was modeled by adding to the existing laminar 
propagation operator, F^, an additional attenuation contribution, 
r c , due to turbulence. 


r t = r L + r 0 = {r Lr + r e ) + ir Ll (17) 

The form of turbulence induced increment in attenuation was derived 
from the form of the attenuation factor of a constant I-R-C, lumped 
parameter, representation of a line. That is, 

T 0 = sifc? (R/Is + if L. (18) 
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V 


Figure 3. Predictions of Attenuation Factor for Turbulent Flow ( Reference 7) 
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V 


Figure 4. Predicted Phase Velocity with Turbulent 

Flow ( Reference 7 ) 


n 


By analogy. 


r c = Re 



(19) 


-wher-e- 


Rt 


2vK N R n 


(20) 


The constants, K and n, were determined by surface fitting the 
turbulent attenuation increment in Figure 3 as a function of Reynold's 
number and frequency. The derived values of these constants are: 


K = 0. 0055 
n = 0. 85 
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III. DISTRIBUTED COMPLIANCES 


III. 1 Distributed Bubbles 


In most practical applications, the liquid rocket propellant in 
the feedline will contain quantities of either dissolved ullage gas or 
propellant vapor, or possibly both. It is assumed that the gases or 
vapor are in the form of small bubbles, homogeneously distributed 
throughout the liquid. The net effect of distributed bubbles in the 
liquid is to modify (reduce) the effective fluid sonic velocity and 
possibly to introduce added damping because of irreversible com- 
pression and expansion of the bubbles caused by a periodic distur- 
bance. 


Two types of entrained gases have been considered in the 
model for the acoustic velocity in a two-phase flow: 

(1) The flow of a single- component, two-phase mixture, 
such as LOX and oxygen vapor, and 

(2) The flow of a two- component, two-phase mixture, 
such as LOX with entrained helium ullage gas. 

For a single- component, two-phase mixture, the actual flow 
process is theoretically bounded by the equilibrium process and the 
constant quality process. In the case of an equilibrium process, the 
vaporization and condensation rates are large enough to ensure that 
the vapor temperature and liquid temperature are always equal; 
both phases are saturated, and the quality (vapor fraction by weight) 
varies only with pressure. Based on the assumption of thermodynamic 
equilibrium, Gouse and Brown® have shown that the speed of sound in 
a single -component, two-phase flow can be expressed. as 

«•- (^r +x ) (s *- s ') KtDJ 5 (2i » 

where X is the mixture quality and (S g - S f ) is the change in entropy 
between the saturated vapor and the saturated liquid state. Evalua- 
tion of Equation (21) requires the extensive use of a temperature- 
entropy chart or a tabular representation of the thermodynamic pro- 
perties of the fluid being considered. 

On the other hand, the condensation and vaporization rates in 
a constant quality process are assumed to be so low that no significant 
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phase change occurs. In a two- component mixture, a constant qual- 
ity exists at -all times if there is no liquid vapor present. Rearrang- 
ing the derivation by Hsieh and Plesset^, the speed of sound in a 
two- component mixture is 


J__ (( p g \ 1 1_ / \ 1 1 1 

c? ” p Up 6 + §p a ) 9a cj e 2 + Vp 6 + ^p^/Ps c 6 3 / 


(22) 


where 


p = mixture density 

P g = gas density 

p a - liquid density 

ca = speed of sound in the liquid 

c g = speed of sound in the gas 

and $ = mass gas /mass liquid. 

As the mass ratio, i, approaches zero, the speed of sound ap- 
proaches that for the pure liquid. Plesset and Hsieh-^ have also 
shown that, for very small gas bubbles dissolved in a liquid, the 
adiabatic speed of sound in the gas should be replaced by the iso- 
thermal speed of sound. The reason is that the gas bubbles have a 
non-uniform temperature distribution in their interior due to the 
finite value of the heat conduction rate. Figure 5 presents the iso- 
thermal versus the adiabatic model for dissolved helium ullage gas 
in LOX. For comparison, the equilibrium model for LOX with 0 2 
vapor is also shown. Reliable experimental data are not available 
to determine which flow process actually occurs for the flow of a 
single- component, two-phase mixture. As can be seen by Figure 5, 
there is a great difference in the acoustic velocity between the equi- 
librium and the constant quality flow processes. These differences 
could lead to a large error in the predicted feedline system re- 
sponse, and more work should be devoted to a study of the actual 
flow process. In the present study, the constant quality model is 
applied to both single- component and two-component flows. A more 
detailed discussion of the acoustic velocity in two-phase flow is pre- 
sented in Appendix B. 

The wave propagation characteristics of the mixture depend 
not only on the gas-to-liquid mass ratio but also on the bubble size 
and disturbance frequency. Bubbles of a given size resonate when 
excited at their natural frequency; this natural or resonant frequency 
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Mass Ratio , <£> 

Figure 5. Variation Of Speed Of Sound With Mass Ratio 
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increases in inverse proportion to the bubble radius* *. An investi- 
gation has been made to determine the variation in effective sonic 
velocity in a liquid with uniformly distributed gas bubbles as the 
excitation frequency approaches the bubble resonant frequency. 


At frequencies approaching resonance. Equation (22) is no 
longer useful, even as an approximation. Spitzer*^ derived an ex- 
pression for the phase velocity at all frequencies by assuming that 
the gas concentration is small. If surface tension is neglected and 
if the bubble diameter is assumed to be larger than 0.001 ft, the 
bubble resonant frequency can be approximated by*^» 


1 / 3T|P 0 _ 44/ \ 8 

° 2ttR 0 \ p i pf Rc f ) 


(23) 


where r\, the poly tropic exponent, is dependent primarily on the 
bubble size and gas thermal conductivity. For very small bubbles, 
such as those assumed to be distributed in the propellant feedline, 
the polytropic exponent is equal to unity. 


The theoretical effect of frequency and bubble size on the 
phase velocity in a mixture of air bubbles in water has been pro- 
grammed and is presented in Figure 6. The mass ratio, §, was 
assumed to be 3.0 X 10 . Note that the speed of sound is only af- 

fected at or near the bubble resonant frequency, which increases 
with decreasing bubble size. For the homogeneous distribution of 
very small bubbles entrained in the liquid propellant, the net effect 
over the frequency range of interest is only to lower the phase velo- 
city of the pure liquid, to the value predicted by Equation (22). The 
effect of larger bubbles, which are assumed to exist only locally 
(cavitation zones), is discussed in Section IV. 


III. 2 Distributed Wall Compliance 

In general, the flexibility of the feedline wall will produce two 
possible effects which may be of concern. First, the wall compli- 
ance will reduce the phase velocity and increase the spatial attenua- 
tion for the longitudinal fluid wave propagation mode (called the 
zeroth mode in Reference 5), and second, the axial wall stiffness 
effect, will permit real wave propagation of the higher order modes 
at low frequency by virtue of a coupling between the fluid and axial 
wall modes. For typical feedline problems, we have concluded that 
the amount of energy which is fed into the higher order modes at the 
line terminations is so small that it can be neglected. Therefore, 
for all practical purposes, the effect of the axial wall stiffness on 
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Figure 6. Effect Of Frequency And Bubble Size 

Water Mixture 






the wave propagation can be neglected. This topic is thoroughly- 
discussed in Section HI. 3, 

The influence of the wall compliance (and mass) on the pro- 
pagation of the zeroth mode cannot-be -neglected-. Gerlach^' investi- 
gated the case of a viscous fluid flowing through a line having elastic 
flexible walls, i. e. , a line whose walls have finite radial compli- 
ance and impedance but do not propagate disturbances axially. By 
equating the radial impedance of the wall to the radial impedance of 
the fluid, the latter being subject to the no- slip condition on axial 
fluid velocity at the wall, the following characteristic equation for 
the eigenvalue problem was obtained: 

(y 3 - g 3 ) p 0 cf 

f s Jipr,) Y £ Ji(kr 0 n = p t hs s + hE t /r 3 (24) 

L P J o (0r o ) ~ k J 0 (kr 0 )J 

The parameters, j3 and k, are separation constants of the solution 
of the Navier- Stokes equations and are related by 

Y® = k 2 + s/v = 0 2 + s 2 /c 0 (25) 

where 

h = tube wall thickness 

r 0 = tube radius 

Et = Young's modulus for tube material 

Pt = density of tube wall 

s = Haplace variable 

c 0 = isentropic speed of sound in the fluid 

V = fluid kinematic viscosity 

p 0 = fluid density 

J 0 , Jj = zero and first-order Bessel functions of the 
first kind. 

Equations (24) and (25) were solved numerically on the computer. 
Figures 7 and 8 illustrate the effect of an "elastic, flexible wall" 

(a wall with radial compliance and mass) on the spatial attenuation 
and phase velocity for the zeroth mode. The area of interest in 
these figures corresponds to dimensionless frequency numbers, 
wr 0 /c 0 , less than 1.0, since it is known that coupled structural- 
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feedsystem instabilities occur in this frequency range. In this range, 
the wall behaves in a spring-like manner and inertia effects are neg- 
ligible; radial attenuation for the elastic wall is virtually indistin- 
quishable from the attenuation of the rigid wall. It has been estab- 
lished that for any practical feedline problem, this is the case; 
hence, only the wall radial compliance need be considered. Wall 
mass effects are negligible. Therefore, the classic Korteweg cor- 
rection to the phase velocity is valid. That is, 

C = / 2P„°co 3 r. \i <26) 

l 1 + Eth J 

It must be remembered, however, that this correction is the result 
of applying a boundary condition to the fluid dynamic formulation. 

III. 3 Axial Wall Stiffness 


In the preceding section, it was stated that the effect of axial 
wall stiffness on higher mode wave propagation could be neglected. 
That statement can be justified by proving the following points: 

(1) Due to observable flow processes in feedlines, there is 

a minimal amount of energy available in the higher modes 
of the fluid for coupling with the wall. 

(2) The energy that is avai lable in each higher fluid mode 
has an attenuation rate, below the cut-off frequency, of 
approximately two orders of magnitude larger than that 
for the zeroth mode. Therefore, the energy level and 
attenuation rate do not permit sustained coupling of the 
fluid and the wall in the higher modes. 

The physical effect of higher modes must be considered, because the 
proposed distributed parameter line model is based on only the 
zeroth mode of propagation, and the existence of substantially 
higher modes would invalidate its functional form. 

Recall that the magnitude of Reynold's numbers in a typical 
feedline indicates a turbulent flow condition. Hence, the velocity 
profile should be relatively flat. Gerlach^ has calculated the pertur- 
bation velocity profiles of the first four modes for the laminar, 
viscous pipe-flow problem. Characteristic results are shown in 
Figures 9 and 10 where F xn (r) is the velocity distribution function. 
Observe closely that, for the wide range of frequencies and damping 
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Figure 9. Axial Velocity Profile Function, F xn (r), for Four 

Modes (cur 0 /c 0 = 1.0, Wr 0 c 0 = 0.01 ) 


22 






\ 

' 


REAL [F xn ] 

1 

-1 

J J 

+1 -1 

0th MODE 1st MODE 

, < 

/ +>1 

1 IMAGINARY 

[ p xn] 

^ 



3658 


Figure 10. Axial Velocity Profile Function , F xn (r), for Four 

Modes (cur 0 /c 0 = 0.2, Wr 0 c 0 s 0.001 ) 
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numbers considered, the zeroth mode velocity distribution is near- 
ly constant across the cross-section. The obvious conclusion is 
that the zeroth mode perturbation profile approximates a turbulent 
profile, and, hence, a minimal contribution from the higher modes 
is required to fill out the turbulent profile. Therefore, it follows 
that the energy level in the higher modes is negligible. 

Based on the foregoing discussion, there is admittedly a 
small contribution of the higher modes to the velocity profile. How- 
ever, Gerlach^ has also shown that the spatial attenuation (the real 
part of the propagation operator) for these higher modes is substan- 
tially larger than for the zeroth mode. The net effect is that the 
higher modes are rapidly damped out over a broad frequency range, 
and therefore have neither the time nor the energy levels required 
to excite (couple) corresponding modes of transmission in the con- 
duit wall. Figure 11 illustrates the elevated damping level of higher 
modes, while the zeroth mode attenuation is roughly insensitive to 
the frequency range. 

The above discussion substantiates the validity of the zeroth 
mode transfer equations for the line. In addition, the foregoing con- 
clusions which were derived on a qualitative basis, have also been 
verified quantitatively by Lin and Morgan^. 
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IV. LOCAL COMPLIANCES 


IV. 1 Cavitation Zones or Large Gas Bubbles 

At elbows or bends in the feedline, or at the pump inlet, local 
pressure drops are likely to create regions of cavitation. In most 
cases, these regions are composed of localized bubbles of vapor. 
These larger cavitation bubbles in the line are modeled by consider- 
ing the bubbles to be local compliances. As the excitation frequency 
approaches the bubble resonant frequency, which can be quite low 
for large bubbles, the volume pulsations may be described by the 
following linear, second-order differential equation: 


Pi 

4rrR 0 


v + bv + 


T1P 0 

Vo 


V = -p(t) 


(27) 


or 


m. 2 V + bv -f kv = -p(t) 


(28) 


where 

Pj g, = liquid density 

P 0 = steady state line pressure 

V 0 = steady state bubble volume 

R 0 = steady state bubble radius 

b = coefficient accounting for viscous, thermal and 
radiation damping. 

The continuity equation relating the flow upstream and down- 
stream of the bubble is 


q 3 - qj. = dv/dt ( 29 ) 

where and q 3 represent the upstream and downstream volu- 
metric flow rate, and dv/dt is the change in bubble volume with 
time. Equations (28) and (29) can be combined conveniently in the 
Laplace transform domain to yield an expression for the change in 
flow rate due to the presence of the bubble. 


Qa(s) - Qi(s) 


- s P(s) 

(m 2 s 3 + b s 4- k) 


(30) 
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where the term l/(m 2 s 3 + bs + k) can he visualized as a '•complex' 1 
compliance. The damping term, b, has been estimated by Devin^, 
and is represented by a combination of the losses originating from 
three processes: 

(1) Thermal damping resulting from the thermal conduction 
between the gas in the bubble and the surrounding liquid; 

(2) Sound radiation damping caused by energy dispersed by 
radiating spherical sound waves when the bubble is ex- 
cited into volume pulsations; 

(3) Viscous damping from viscous forces at the liquid-gas 
interface. 


In the literature, there is considerable disagreement as to 
the theoretical expressions that should be used to evaluate the 
damping constant, 6, which is related to the coefficient of damp- 
ing, b, by 



(31) 


The reason for this disagreement is that, while the theory is appli- 
cable to all frequencies, experimental verification has been obtained 
only at the resonant condition. However, it is hypothesized that the 
theory which is summarized here and has been incorporated into the 
computer code is valid over a broad range of frequencies. 


The total damping constant is given by 

5 = 6th + ^rad + ^vis 
or, b = bth H* brad + bvi B 
where 

/ sinh (20i R 0 ) + sin (20!R O ) 1 > 

_ I cosh(20iR o ) - cos(20jR o ) 0iR o ( RPq 
th ~ ) 20x R 0 sinh (20 ! Rq) - ain(20 1 R o ) ( V 0 w 
V3(Y-1)* cosh(20!R o ) - cos(20iR o ) ) 


brad 


PjUJ 3 

4ttc^ 



(32) 

(33) 

(34) 

(35) 
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In these expressions. 


(i _= liquid vis.cosity 

■Pj i = liquid density 

Cj ^ = liquid speed of sound 

Y = ratio of specific heats irrespective of the' thermo- 
dynamic process indicated by the polytropic expo- 
nent, T). 


Also, the argument 


0iRo 



Ro 


(36) 


is a measure of the rate at which heat is conducted over a distance, 
R 0 , from the bubble center to the liquid-gas interface. The thermal 
diffusivity of the entrained gas is denoted by Eh . 


For zero damping of the oscillations, 

Q x (s) - Qg (s) =- , sp ( s > . 

k *+?«•) 


(37) 


At frequencies considerably below the bubble natural frequency, U) 0 , 

s P(s) 


Qi(s) - Qa (s) = 


(38) 


where 0J o is defined as 

U) 0 = 'sfk/ra^' (39) 

For these low frequency cases, the compliance becomes simply 
Vo 


c - 


11 Pc 


(40) 


The correct value of the polytropic exponent r\ in the stiff- 
ness term k = 1) P 0 /V 0 depends on the bubble size and frequency. 
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as well as the thermal properties of the gas and liquid. This varia- 
tion with bubble size and frequency of excitation is shown, in Figure 
12. For a given bubble undergoing successive compressions, the 
gas near the center behaves adiabatically, while the gas near the 
liquid-gas interface undergoes no change in temperature since the 
liquid behaves as a large heat sink. The value of the polytropic 
exponent has been found to be 


11 = 


Y [ 1 + &th 3 ] " X 

r 3 ( Y- 1 ) f sinh (20! R 0 ) - sin (20iR o ) \ “I 

L * 20 1 R O \cosh(20 1 R o ) - cos(20!R o ) / J 


(41) 


where 6 tb is equal to the term in brackets in Equation (33). 


IV. 2 Complex Side Branch 

The types of side elements typical of a rocket feed system 
(accumulators and gas-filled pressure sensing lines) behave as local 
compliance -like elements. 


In general, any type of complicated side element may be 
modeled; these side elements can be visualized as local or lumped 
"compliances, " even though they may not be compliance^-like ele- 
ments. The model used for an arbitrary complex side element, 
such as a gauge and line connected onto a feedline, may be put in the 
form of a local compliance, as shown in Figure 13. This element 
is assumed to consist of a side branch line having an inertance, I, 
and resistance, R, and with a capacitance termination, C. The 
line inertance, I, and the line resistance, R, have been pro- 
grammed respectively as: 


I — P o R/ A b 


and 


R = 


128 


1-tR 

TTdb 4 


(42) 


(43) 


Here, p 0 = fluid density, L = side branch length, and A b 
is the branch cross sectional area; p is the fluid viscosity, and d b 
is the side branch diameter. To consider the effects of compressi- 
bility in the side branch, the inertia and resistance effects are neg- 
lected, and 
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POLYTROPIC EXPONENT 



Figure 12. Polytropic Exponent for 0 2 Bubble in LOX vs Angular Frequency and Bubble Size 



(a) PHYSICAL SITUATION (b) MODEL 



(c) MODEL EQUATIONS IN LAPLACE DOMAIN 

CONTINUITY IN LINE 

Q,(s ) - Q 2 (s) = Q 3 (s) 

GAUGE LINE ( Assumed simple inertance plus resistance only ) 
Q 3 (s) «Q 4 (s) 

P 2 (s)- P 4 (s ) = IsQ 3 + RQ 3 
CAPACITANCE 

Q 4 (s ) = C$ P 4 
"COMPLIANCE MODEL" 

Q,(s) - Q 2 (s) ic s 2 + RC S + r} S 

Figure 13. Example Treatment Of More General 
"Local Compliance" in Feed Line 
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* = *p (It) for gases 

q* = ~ (|jr) for liquids (44) 

where qj. is the flow into the branch, V is the volume of the com- 
pressible fluid contained within the side branch, Y is the ratio of 
specific heats for the fluid (if a gas), and H is the liquid bulk 
modulus. 

For the case of a gauge connected to a feedline containing a 
liquid by a branch line partially filled with air, the lumped model is 
obtained by summing the effects of the branch inertance, resistance, 
and capacitance. 

The resultant expressions in the Laplace domain relating the 
pressure and flow upstream and downstream of the side branch 
become 


Fs(s) = Pi(s) 


(45) 


Q 2 (s) - Ql(s) - [ ICs S + RCs + l J sP l( S ) (46) 

The quantity in the brackets can be visualized as a complex 
"compliance"' for this case, and represents a local compliance 
so far as the feedline is concerned. 
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V. COUPLED RESPONSES 


Four possible modes of structural-hydraulic coupling are 

shown in Figure 14. In cases a, c, and d, an externally imposed, 

rigid-body acceleration acts on a length of feedline. Because of 

viscous shearing forces at the propellant- conduit interface and/or 

body forces applied through the end impedances, the imposed line 

motion will couple with the fluid motion to produce modifications to 

the pressure-flow equations for the stationary line (Section II. 1). 

The functional forms of these modifications are derived herein. 

\ 

V. 1 Accelerated Line (Cases c and d) 


The fluid may be 
linearized, perturbation 


characterized by the familiar axi symmetric 
equations of motion, continuity and state. 


dv 

at 


1 _ 3p rjfv 1 3v ] 

Po dx + V Ld r 2 + r dr J 



(47) 


dp 

at 


+ 



0 


(48) 


dp _ dp 

Po = K ’ 


K — P o C o 


(49) 


Implicit in this formulation is the assumption that the radial pertur- 
bation velocity is negligible compared to the axial perturbation velo- 
city, v. Note that the pressure has not been assumed to be inde- 
pendent of the radial coordinate. The body force F x , can be re- 
placed by the equivalent D'Alembert force, ~p 0 a x , due to the im- 
posed acceleration. Utilizing this substitution, combining the con- 
tinuity and state equations and applying the Laplace transformation 
with zero initial conditions, yields the following pair of coupled 
differential equations: 


sV 



dp 

dx 


+ V 


rd 3 v 
L dr 2 




(50) 


PoCp 2 dv 

s dx 


(51) 


where P and V represent the transformed fluid pressure and velo 
city. These equations can be combined into a single equation in the 
dependent variable V. 
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Figure 14. Elementary Structural-Fluid Coupling Cases 
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2 


(52) 


0_ 

S 


d 2 V f^V 

— r + v ^3" + 
3x 2 L Sr 3 


1 _ 

r 



Making use of the fact that A x is a function only of the Laplace 
variable, s, and redefining the dependent variable to be 


U (r, x, s) = Y + A*/ s 


(53) 


results in the following easily solved. 


ordinary differential equation: 


v r a 2 u I TT "| , 5 2 u 

s LSr 3 + r 5r v J s 2 bx 2 


0 


(54) 


This equation is amenable to solution by the standard separation of 
variables technique which yields 

U = (A cosh Xx + B sinh Xx) J 0 (ar) (55) 


where 



X 2 - separation constant 
A, B = integration constants. 

By definition, X is the propagation operator, Y(s) which governs the 
spatial attenuation of axially propagated waves in the fluid. This 
operator may be stated a priori: 

V(S) ' 77 f 2jl(Sr„) \j < 56 > 

\ ?r 0 J 0 (?r 0 ) J 

where 

= -s/v , r 0 = conduit radius 

In. addition, the characteristic impedance of the line, Z 0 , is defined 
to be 


Z c = 


P o Co 


S 


Y(s) 


(57) 



Utilizing these definitions, the fluid pressure and velocity become 


V(r,x, s) = [A cosh Yx + B sinh Yx] J 0 (Ctr) - A x /s (58) 

P(r,.x,.s) .= -Z 0 .[A sinh Yx + B cosh Y-x-] J- 0 -(ar)- (59)' " 

The dependence on the radial coordinate may be eliminated by 
averaging these expressions across the cross section, 

— 2 

V = [A cosh Yx + B sinh Yx] Jj. (ar 0 ) - A x /s (60) 

— 2 

F=~Z 0 [A sinh Yx + B cosh Y x ] Ji (ar 0 ) (61) 

ocr 0 

Constants, A and B, can be evaluated at the origin of the line (z=0) 
by applying the boundary conditions. 

V (o, s) = Vi ; P (o, s) = Pi (62) 

which results in the following expressions: 

V = Vx cosh Yx sinh yx - A x /s (1-cosh yx) (63) 

Z c 

P = Pj cosh Yx - ZeVi sinh yx - — — K ' sinh Yx (64) 

s 

The desired four-terminal representation for the mean exit pressure 
and flow velocity are obtained by evaluating the above expressions at 
x = L. 


V 3 = Vi cosh yL - Pi/Z c sinh yL - V x (1-cosh yLi) (65) 

P 2 = Pi cosh yL, - Z c Vi sinh yL - V x Z c sinh YL, (66) 

where the acceleration has been replaced by its Laplace equivalent 
s V x . These equations can be referred to local volumetric flow rate 
by introducing the line cross-sectional area, A. 

V. 2 Relative Motion of Bellows and PVC Joints 

Bellows 

Bellows are commonly used elements in propulsion feed sys- 
tems, and in their most important application, they serve as a "fix' 1 
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for the POGO instability. Bellows vary in their structural com- 
plexity from the simple array of axisymmetric convolutes shown in 
Figure 15a to the configuration in Figure 15b which incorporates a 
gas trap liner. The following discussion is directed toward the 
latter configuration which includes the simple bellows as a special 
or degenerate case. 

When the ends of a bellows execute relative motion, the change 
in volume flow rate between the inlet and exit can be considered to be 
the sum of two separate effects: 

(1) A flow rate change due to the compliance of the gas 
trapped under the liner, and 

(2) An apparent local volume production. 

It should be noted that the first item exists even in the ab- 
sence of relative motion. 

First, the pressure drop across a bellows in the flow direc- 
tion is expressed in the Laplace domain by 

P 2 (s) = F(pV|,G) Pi (s) , (67) 

where F(pV|, G) is a loss factor, dependent on the dynamic pres- 
sure and bellows geometry. Now, the change in flow rate due to 
trapped gas can be written in terms of the bellows inlet pressure, 

qs' - qi 7 = - C 

or in the transform domain 

Qa'(s) - Qj/{s)=-Cs Pi (s) (69) 

where C is the lumped compliance of the gas and is defined by the 
ratio of the gas volume to the gas bulk modulus of elasticity. 

The apparent volume production can be derived from the con- 
tinuity equation 

q 3 " - q^' = dv/dt (70) 

or in the Laplace domain 

Qa'(s) - Qi"(s) = -sVol(s) (71) 
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a. Simple Bellows 



( )' - quantities associated with compliance 
of trapped gas 


( )" - quantities associated with volume 
production 

b. Bellows With Gas Trap Liner 

3442 


Figure 15. Nomenclature For Relative Motion Of Bellows 
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The volumetric production is related to the axial displacement of 
the bellows ends by 


Vol(s) = Kb [Xj, (s) - Xa(s)] (72) 

where Kb is a "volume change" constant, also dependent on the 
bellows geometry. The resulting flow equation is 

0/(3) - Qi"(s) = sK b fXx(s) - X 2 (s)] (73) 

or in terms of relative velocities 

Q/(s) - Q x "(s) = K b [Vi(s) - V s (s)] . (74) 

The combined effect, expressed as the sum of Equations (67) and (74) 
produces the following result for the change in flow rate: 

Qa(s) =.Qi(s) - CsPi (s) + Kb [V t (a) - V s (s)], (75) 


Equations (67) and (75) describe completely the four terminal pres- 
sure-flow relationship for a bellows. 

PVC Joints 


In many applications, the PVC (pres sure -volume compensa- 
tor) joint is often used in place of an elementary bellows where the 
local volumetric change effect described above is judged undesirable. 
A typical PVC joint, shown in Figure 16 , generally consists of a 
bellows to allow line flexibility (either axially, or to allow angula- 
tion, or both) and a compensation chamber to "absorb" . the excess 
volume of fluid created in the line when the bellows is displaced 
axially. In theory, the PVC joint is designed to permit axial length 
changes without introducing a corresponding fluid volume disturbance 
in the line. In practice, however, the PVC joint is not a perfect 
compensation device. 

The "volume flow production" of the bellows Q 4 , may be 
written as 

Q 4 (s) = Qg (s) - Q 3 (s) = K b [Vj. (s) - Vg (s)] (76) 

where Vi(s)*, i= 1 or 2, is the transformed velocity associated with 
coordinate Xi . Similarly, the ideal volume compensation of the PVC 
may be written as 
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Figure 16. Illustration Of Typical PVC Joint 
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Q s {s) = Q t (s) - Q 3 (s) = K c [Vj. (s) - V 2 (s)]. 


(77) 


Continuity requires that 

Qs(s) = Qi(s) + Q 4 (s) - Q 2 (s). (78) 

Combining Equations (76), (77), and (78), the resultant flow be- 
comes 


Qs(s) = Qi(s) + (K b -K 0 ) [Vj, (s) - Yg (s)] . (79) 

It is normally intended that K b = K 0 so that perfect compen- 
sation is obtained. However, there is a dynamic effect associated 
with the action of the compensation unit, so a more general form of 
Equation (77) is 

Qs (s) = ~~ [Vi (s) - V 3 (s)] (80) 

where G(s) might be a simple lag, e» g. , (1 + ts), caused by a com- 
bination of the PVC structure elasticity and fluid. compressibility. 

Equation (79) was selected as the PVC flow rate expression in the 
feedline computer code. The pressure drop across the joint was 
assumed to have the same functional form as Equation (67). 

V. 3, Forced Changes in ‘Line Length 

The POGO phenomenon is a coupled dynamic instability involving 
the vehicle structure, propulsion system, and propellant feed system. 
Propellant feedlines are normally supported by the vehicle structure at 
discrete points on the line as opposed to a continuous support. A forced 
change in the length of the line between consecutive supports occurs when 
vehicle structural inputs to these supports differ in phase and/or magni- 
tude. In this section, the effects of changes in line length are analyzed 
and the transfer matrices that describe the pressure flow response to 
this type of excitation are derived. 

The effect that length changes produce on the pressure-flow re- 
lationships for the liquid in the line can be analyzed from two different 
viewpoints. The first approach requires the solution of the fluid dynamic 
equations of motion subject to an inhomogeneous boundary condition on the 
spatial variation of the axial wall velocity as dictated by the structural 
inputs, Vj, and V 2 . An alternate approach assumes that the dynamic 
length change can be modeled as the linear superposition of (1) the axial 
vibration of the rigid line, and (2) a volume production region that reflects 
the relative motion between the ends of the line as indicated in Figure 17. 
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Figure 17. Illustration Of Approximate Modeling Procedure For 

Forced Changes In Line Length 
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The first method of analysis begins with the Laplace transform 
version of the linearized (1) radial and axial equations of motion for 
the fluid, (2) continuity equation, and (3) the equation of state. These 
equations are reduced to one -dimensional form by averaging across the 
feedline cross section. The solution of the resulting differential equa- 
tions for the transformed axial velocity is matched to the transformed, 
inhomogeneous compatibility condition on wall velocity. This boundary 
condition is obtained by solving the undamped wave equation for the 
wall. 


The physical problem is shown in Figure 18. ' Liquid pro- 
pellant is flowing through a constant area duct of length L, internal 
radius r 0 , and wall thickness t. The terminal ends of the line are 
acted upon by structural velocities having arbitrary relative phase 
and magnitude. The relative motion due to this excitation is assumed 
to be small in the linearized sense. 


1 

— V, (t) 


r, 

J. SSSSSS/SS/ss.ys //// *7* 


> — ^ 

r o 


9 — A 


7'/;/;/// / / s' / / / / ZZZL 


x=o x = L 


V 2 (t) 


p=p <t) p=p 2 (t> 

q=q, (t) q =q 2 <t) 


Figure 18. Geometry For Forced Changes In Line Length 

The following assumptions and assertions are made with 
regard to the fluid flow: 

(1) The flow is viscous, incompressible and axi symmetric; 

(2) Viscosity remains constant; 

(3) The flow is laminar; this assumption is not a severe 
constraint because the added attenuation due to turbu- 
lence can be adequately taken into account through the 
use of the modified propagation operator described in 
Section II; 
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(4) The pressure and the axial and radial components of 
flow consist of a steady- state level plus a perturbation 
value; 

(5) The axial velocity perturbation, ,v.,. is much larger 
than the radial perturbation, v r ; 

(6) Velocity gradients in the axial direction contribute a 
negligible amount to the viscous dissipation forces; 

(7) Spatial derivatives of the density can be neglected since 
the compressibility of the liquid propellant is extremely 
low. 


D'Souza and Oldenberger*^ have shown that when these assumptions 
are applied to the Navier- Stokes equations in cylindrical coordi- 
nates, the continuity equation and the equation of state, the follow- 
ing linearized differential equations describe the flow field: 


dv dp rd s v 1 3v 1 

Po d7 = " dx + ^ L d^ + 7 d7 J 


(81) 


1 dp dv r v r dv 

— + - + = 0 
K dt dr r dx 


(82) 


These expressions are easily recognizable as the equation of motion 
in the axial direction and the combined continuity and state equations, 
where n. is the fluid bulk modulus of elasticity. In addition, ab- 
sence of an equation of motion in the radial direction implies that 
the pressure is constant across the cross section, i. e, , p = p(x, t). 


The relative motion in the line wall is governed by the wave 
equation 


d 2 y _ 2 5 S y 
dt 2 " Ctf dx 2 ’ 



(83) 


where y is the axial displacement of a point in the wall and c w is 
the longitudinal wave speed in the wall. Equation (83) is subject to 
the imposed boundary conditions 

^ = and = (84, 

The no- slip compatibility condition at the wall-fluid interface is 
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(85) 


dy(x, t). 

at 


v(x, r 0 , t) 


The solution for the flow field proceeds in the Laplace trans- 
form domain noting that the average pressure perturbation at an 
arbitrary cross section, p(x, t), is equal to the actual perturbation, 
p(x, t). The transform of the equation of motion with zero initial 
conditions yields * 


a 2 v i av s / j_ ap\ 
dr 2 + r Sr v \ + PoS dx J 


0 


( 86 ) 


where the upper case letters indicate transformed variables. Since 
P is independent of the radial coordinate. Equation (86) can be re- 
duced to an expression involving only one dependent variable by the 
transformation of variables 


U = 



dp 

3x 


(87) 


Combining Equations (86) and (87) yields a Bessel type differential 
equation with one finite solution at r = 0. 

U = f(x, s) J 0 (?r) (88) 

where §= \/-s/v and f(x, s) is an undefined function which must be 
determined. Therefore, the transformed velocity, V, becomes 

V = f(x, s) J 0 (Sr) - (8.9) 


Solution of the wave equation in the wall subject to the prescribed 
boundary conditions yields 


Y(x, s) = Bi 


sx 


cosh — -f B 3 sinh 


sx 

c w 


where B* 


Vi(s) 


and B s 


V 3 (s) - Vi(s) cosh sLi/cy 
s sinh sL/ c w 


(90) 


Combining the transformed compatibility condition. Equation (85), with 
(89) and (90), produces 


1 dp 

p 0 s dx 


, sx „ . , sx 1 

cosh • — • + B 3 sinh — 
C* Cy J 


(91) 
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Substituting this result into Equation (89) and averaging over the 
cross section yields 

V= f(x, s) p J ^ -— - J 0 (§r 0 )1 + s[bi cosh— -t B 2 sinh— 1 (92) 


As before, barred quantities indicate average values. If the con- 
tinuity equation is transformed and then averaged across the cross 
section, the result is 


p = .1 

s 3x 


(93) 


Note that the dependence on radial velocity vanishes since there is 
no net radial flow. Differentiating this result with respect to x and 
incorporating Equations (92) and (91) yields a second-order, in- 
homogeneous, ordinary differential equation for the arbitrary func- 
tion f(x, s). After considerable algebraic manipulation, the equa- 
tion takes the following form: 


d 3 f 

dx? 


-Y 2 f 



sY 2 

J 0 (Sr 0 ) 


g(x, s) 


(94) 


where 


Y 2 = 



2Ji(gr o) 1 
§r 0 J 0 (§r 0 )J 


and 


g(x, s) = B x cosh 


— + B 2 sinh 
c* 


sx 

c w 


A particular integral of Equation (94) is easily obtained by the 
method of undetermined coefficients because of the cyclic nature of 
the inhomogeneous function g(x, s). The general solution of Equa- 
tion (94) is 


f(x, s) = (Ax coshYx+A 2 sinhYx)+ [Q. cosh— + Qgsinh-^ ) (95) 

V c *r 


:inh— ) 


where Ax and A 2 are undetermined functions of the Laplace variable 
s and Cx and C 2 are related to Bx and B 2 by 


Ci = 


-V s )jo(?r.) 


i = 1,2 


(96) 
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After simplification, the transformed velocity field. Equation (92) 
becomes 


V= (Aj, coshYx+ A 3 sinhYx) J 0 (?r 0 A + 



Noting that k = Pt>c 0 S , the pressure field is obtained from Equation 
(93). 



Coefficients A x and A 2 can be eliminated by evaluating the pressure 
and velocity at the extremities of the line, i. e, , 


P{o, s, } = P x P(L, s) = P 3 

( 99 ) 

V(o, s, ) = V x V(L, s) = V 3 

In the process of making these evaluations, the characteristic impe- 
dance, Z c , can be defined in terms of the propagation operator, Y, 


7 PqCq S Y 
c “ 

s 

Applying boundary conditions (99) to Equations (98) and (97) 
incorporating the expressions for B x and B 2 , and imposing the 
condition that V 3 (s) = G(s) Va.(s) and Q(s) = AV(s) 


r 

0* 




cosh YE - Z„ sinhYL 



(-4 - v 3 ) 

"/ C S 3 (^(S) 

~v 

l U 2_ 


_-Z 0 1 sinhYE coshYL 



(tH 

_<* 2 _ 


where 
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Oti = 


P O Cp 


ft sinh^. 

\ C W 


Ysinh yL ) + 


v 


(G-cosh—^ 


sinh 


sL 


— Xlcosh^-^- cosHyL) 
c H \ c H y 


(102) 


Os = - 


^cosh - ^ 1 - coshyL] 

i C H 




[G-cosh— j, 

V 2*1 (r* 


sL 

C w 


sinh-^--^— sinhyLi 
c w c w Y 


') 


(103) 


The assumed dependence of V 2 (s) on Vj^s) is entirely realistic since 
the vehicle structural velocities, which are the forcing functions, are 
generally describable by transfer functions. 

Equation (101) represents the desired result for forced changes 
in line length. The column matrix in Equation (101) represents a 
pressure-flow modification to the square matrix that describes pres- 
sure and flow in the absence of external excitation. It can be shown 
that when Vj, = V 3 the above expressions degenerate to 


Ps" 


coshyL -Z C A 1 sinh YL 


Pi 


Z c sinhyL 


— 




+ Vf 


Q 2 


-AZ a sinhYL coshYL 


Qi 


A(l~ coshyL) 


where 


(104) 


£ _ 2 J~i ( ?r 0 ) 

?r 0 J 0 ( §r 0 ) 

This latter result was predicted by Gerlach^ for rigid body, axial 
vibration of a line, where coupling occurs through viscous shear. 


In the second modeling approach, the pressure-flow equations 
for the rigid body motion of a line are modified with a volume produc- 
tion element that reflects the relative motion of the ends of the line. 
For a rigid line oscillating with velocity, Vj. , 


*Ps 


coshYL -A 1 Z 0 sinh YL 


Pi" 


Z c sinh yL 

_ / 

= 

_ -1 . 



-Vi 


Q2 


-Z c A smh YL coshyL 


Qi, 


A(l-cosh yL) 


The rate of volume generation due to relative motion is 
Q/ = A(V a - Vi) 


( 106 ) 



Total volumetric flow rate is then the sum of Q 3 / and Q^. 
bining Equations (105) and (106) yields 


Com- 




Pi Z c sinhyL 

-V x (107) 

Qj. A^2-cosh yL-G(s)j 


As before, G(s) describes the relationship between the externally 
imposed velocities, V x and V 2 . 


Equations (101) and (107) cannot be compared directly since, in 
the first case, the excitation is coupled to the fluid by a boundary 
condition and, in the latter case, the coupling takes the form of a 
body force applied to the fluid. 


We believe that the result obtained by treating the problem 
as an inhomogeneous boundary valued problem is the most accurate, 
and it is this expression. Equation (101), that has been incorporated 
into the computer code. 


V. 4. Mounting Stiffness 


Propellant feedlines are anchored to the primary vehicle 
structure by mounting brackets that exhibit varying degrees of 
elasticity and damping. These brackets act as filters which modify 
the structural excitation that is transmitted to the feedline. There 
are two techniques that can be utilized in modeling the effect of 
mounting stiffness: (1) discrete parameter-acceleration, and (2) 
impedance. Both of these methods are summarized in this section. 
The transfer equations for the discrete parameter-acceleration 
technique have been programmed as Subroutine Eight in the com- 
puter program while the impedance approach is programmed in 
Subroutine Ten. 


Discrete Parameter-Acceleration Technique 


Consider the idealized configuration shown in Figure 19. The 
entire line, which has two bends, is assumed to be infinitely rigid and 
elastically restrained. A single linear spring in parallel with a vis- 
cous damper represents the combined effect of all the discrete 
mounting stiffnesses whose line of action coincides with that of the 
equivalent mounting stiffness. In Section V. 1, an expression was 
derived for the pressure and flow rate in a vertical line that is ex- 
cited by a velocity, V^, applied directly to the line. That expres- 
sion forms the starting point for the analysis of the present configura- 
tion. Initially, the velocity, V^(s), is replaced by its Laplace equivalent, 
ajp(s)/s. It was shown in Section V. 1 that the four-terminal represen- 
tation of a vertically accelerated line is 
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Figure 19. Model For A Line With Mounting Stiffness 



Ps" 

= 

"coshYL -Z c A 1 sinhYL _ 


'Pa" 

a^(s) 

s 

Z 0 sinh YL 

_Qs. 


-AZ c 1 sinhYL coshYL_ 


_Qa_ 

1 

A{1- coshYL) 


( 108 ) 


The problem is to express the resultant line acceleration, a^{s), in 
terms of the applied structural acceleration, ai(s), the viscoelastic 
support parameters, and the fluid mass contained in the horizontal 
limbs of the feedline. By neglecting the elbow at the exit to the fuel 
tank, this aspect of the problem is easily overcome by viewing the 
entire system as a damped, harmonic oscillator whose mass is equiva- 
lent to the combined liquid mass in the horizontal limbs. The support 
is excited by a^s), and the inertial effects in the vertical line seg- 
ment are taken into account by applying the dynamic fluid pressure 
forces on the projected areas of the elbows on the system mass. 

During oscillatory motion, the mass of fluid contained in the 
transverse limbs of the feedline contributes an added mass, (m 1 + m 3 ), 
effect. The inertial loading, due to the fluid mass in the vertical seg- 
ment of the line, is replaced by the equivalent dynamic pressure forces. 
For a support of negligible mass, the differential equation of motion 
for mass M about its equilibrium position is. 

Mx = -k{x-y) - b(x-y) + A(p 3 -p 3 ) (109) 

Transforming this result into the Laplace domain and noting that 

£ {x} = X, ^Cy}=Y (110) 

and 

s 3 X = ajg(s), s 2 Y = ai(s) (111) 

the following result is obtained. 

a^<s).a,(s)[ 5s ^-] t A[P 3 (s)-P a ( s )][ ^ 5a ; s 3 Mk ] (112) 

The first term in Equation (112) reflects the externally imposed 
acceleration, and the second term reflects a passive inertial loading. 
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A straightforward combination of Equations (112) and (108) produces, 
after some manipulation. 


■p 3 ~ 


r Pqi+ft 7 

J=qs . 1 
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where the Eij (i, j = 1,2) are the elements of the original square 
matrix in Equation (108) and 


a 


sb+k 


s(Ms 3 +sb+k) 

As 

(Ms 2 +sb+k) 


Z 0 sinhYL 


Z e sinh YE 


a 


- .(Jw> A(1 - co3hvL| 

A 3 s(l-cosh YL) 

Ms 3 +sb+k 


(114) 


All other terms have been defined previously. The input accelera- 
tion, ai(s), could be replaced by the equivalent velocity represen- 
tation, s Vi(s). However, this replacement would require a minor 
computer code modification since the program currently determines 
the response to acceleration, ai(s). As the support becomes pro- 
gressively stiffer, Equation (113) degenerates into the expression for 
a simple accelerated line. 


Equations (113) and (114) represent the complete transfer 
functions for the discrete parameter -acceleration technique for a 
feedline segment with external accelerations applied through a mount- 
ing structure with arbitrary stiffness. This technique has been pro- 
grammed into the computer code as Subroutine Eight, 

A second method of analysis of problems involving structural 
mounting stiffness will now be presented, where the line accelerations 
(or velocities) are expressed in terms of the structural driving point 
impedance. 
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Impedance Technique 


The impedance technique will now be applied to the same line 
configuration that was used in the previous case with the exception 
that the mounting stiffness which consisted of a grounded acceleration, 
3-i, and a viscoelastic support is replaced by a driving point impe- 
dance, Z g . As before, the starting point for the analysis is the four- 
terminal representation of line that is subject to velocity excitation, 
Vjg(s). It is required that the line velocity be expressed in terms of 
the imposed forces and ultimately in terms of the driving point im- 
pedance. The support exerts a force, F, on the feedline, and the 
equation of motion for this system that corresponds to Equation (109) 
is 


Mxjj = -F + A (pa-pa) (115) 

Applying the Laplace transformation to this expression yields 

Ms V^(s) = -F(s) + A [p 3 (s) - p 3 (s)] ( 116 ) 

The driving point impedance is, by definition, 

EM 


Zs(S) " Vjt(s) 

Eliminating F(s) between Equations (116) and (117) produces 
A(P 3 -P 3 ) 


(117) 


Vje(s) = 


Zg + hi s 


(118) 


Equation (118) may be introduced into Equation (108) and the follow- 
ing result is obtained after rearrangement: 
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AZ cS sinhTg 
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A 3 (l-coshr 3 ) 2 
Zg -f- hi s 


(119) 


Lij = elements of original square matrix in 
Equation (108). 
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The obvious difference between Equations (113) and (119) is 
that the latter does not contain the trailing column vector that re- 
flects the enforced excitation. In the latter formulation the effect 
of the mounting support (structural impedance) is absorbed into the 
square matrix. Both formulations degenerate into the simple sta- 
tionary line representation for an infinitely stiff support. 
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VI. GENERALIZED FEEDLINE COMPUTER. CODE 


A versatile computer code has been generated to calculate 
the frequency response of propellant feedlines. Program philosophy 
and logic are discussed in this section, and the mechanics of setting- 
up and execution including specific example problems are presented 
in Appendix D„ 

The construction of a computer code to generate the frequency 
response of a propellant feedline is a relatively straightforward task 
when the sequence of feedline components {lines, bellows, PVC joints, 
etc. ) is fixed by a specific design configuration. However, for pre- 
liminary design analyses, it is advantageous to determine the effect 
that a reordering of line components would produce on the overall 
feedline response. For example, various sequences of components 
could materially affect the structural-hydraulic coupling that is 
known to influence the instabilities that were delineated in Section I. 

It can readily be visualized that reprogramming the feedline re- 
sponse equations for numerous changes in the component sequence 
is an extremely expensive, time-consuming task. To alleviate this 
situation, a master computer code was written in FORTRAN IV for 
the CDC 6400 to provide the systems engineer with a design tool for 
obtaining the frequency response of a feedline in which the type, 
number and sequence of basic line components between the propellant 
tank and the turbopump inlet can be specified in a completely arbi- 
trary fashion. The dominant criterion employed in generating this 
code was that the "user" be required to perform a minimal amount 
of manual conditioning of a given problem prior to machine execution. 

Formulation of the code is based on the fact that the four- 
terminal pressure-flow relationship across any line component can 
be represented in matrix form when the perturbation equations are 
transformed into the Laplace domain. That is, with no loss in 
generality. 




■Pf 


= Di 


-Qi+i. 


_Qi. 


± KCd i = 1,2,. 


..n (120) 


where n represents the number of components in the line. In this 
expression, the transformed output pressure and flow for a given- 
component are related to the corresponding input quantities through 
a 2X2 square matrix, D j , plus a column matrix, Cj , that is 
present only if that particular component is being excited by an 
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external forcing function, K. The parameter, K, may be the 
Laplace transform of an externally imposed acceleration or velo- 
city, or the volume perturbation due to a side branch pulser. In 
general, it is desired to determine the perturbation pressure at 
some— point -in- the- line -in re spons e- to -any one of- the- -e-lass- of -forcing 
functions that K represents. With regard to propellant feedline 
analyses, the important pressure point is located at the entrance to 
the turbopump because dynamic variations in the inducer inlet suc- 
tion head are instrumental in the growth of structural-fluid dynamic 
instabilities. Having established the desired output variable, the 
functional form of the transfer function for the overall line is ob- 
tained by applying Equation (120) to each line component followed by 
sequential matrix substitution to arrive at the following generalized 
transfer equation: 


P - D Q ±KB (121) 

Through experience we have found that matrix D. can be stated 
a priori for any feedline composed of n ordered components. 

D = D n D n _ x D n _ 2 B, (122) 

Matrix 33, however, does not possess such a well-defined property. 

In general, 

B = B x + B 3 + +B b (123) 


The value of "m" is strongly dependent on the line configuration as 
is the matrix structure of each Bj . The distinguishing feature of 
each Bj is that it consists of a column matrix of one of the C^'s 
pre-multiplied by one or more Ill's. Matrix 3? in Equation (121) 
contains the desired pressure response to the perturbation, K, and 
matrix Q contains the pressure-flow perturbations at the exit to 
the fuel tank. By assuming that the impedance at the tank exit, Zj , 
is known and that the turbopump -injector- engine combination can be 
lumped into an equivalent impedance, Z t , Equation (121) can be ex- 
panded into its constituent equations to obtain the explicit form of 
the transfer function 

Pn ± Ijd gJix+dgsl >11 - ( d U Z 1 + d l3) fr&J, (124) 

K (dg^Zi + dss) - (dijZi + d i2 ) / Zt 

where the and dij are the elements of matrix 13 and D, 

respectively. 
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To illustrate the concepts that have been presented, consider 
the feedline shown in Figure 20. It is desired to determine the 
transfer function relating the pressure at the turbopump inlet, P s ., 
to the line excitation velocity, Vjj. Adjacent to each line compo- 
nent is a matrix expression of the type presented in Equation (120). 
Performing the previously indicated matrix substitution yields 


Ps 



'ZiQT 


= D 5 Eh 

D 3 Dj_ 


Ps/Zt. 



. Qx _ 


V x (Ds C4 D5 D4 C3 - Ds D4 D3 C 3 ) 


(125) 


The subscripts on pressure and flow rate pertain to specific points in 
the line while the subscripts on matrices Di and refer to the 

component location in the line. Equation (125) completely defines the 
generalized matrices in Equation (121). An expression similar to 
Equation (124) follows directly from Equation (125). 

The foregoing discussion forms the basis for constructing a 
computer code to determine the frequency response of a feedline con- 
taining an arbitrary number and sequence of components. The com- 
puter code contains a controller program and a separate subroutine 
for each component in which the elements of each D* and Ci are 
calculated. The program has the capability of synthesizing a feedline 
with as many as fifteen different types of components. However, the 
program listing that is presented in the Appendix contains eleven 
component subroutines; the four vacancies are available for the 
addition of components during future programs. In addition, the 
source deck contains two subroutines for constructing matrices jB 
and _D and one subroutine for calculating the Bessel functions J x 
and J 0 and their ratio for complex arguments. 

Another subroutine calculates the speed of sound in a liquid 
that may or may not contain a homogeneous distribution of bubbles. 

In either case, the radial elasticity effects of the line may be taken 
into account through the Korteweg equation. 

Since substitution of iu> for the Laplace variable "s" implies 
calculations in the complex plane, liberal use has been made of the 
complex variable operations afforded by FORTRAN IV. Complex 
quantities have been subscripted with as many as three indices to 
facilitate bookkeeping within the program. 

In summary, the user merely defines the feedline structure 
as was done in Figure 20 and indexes the components sequentially be- 
ginning with the component that is attached directly to the propellant 
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Figure 20. Matrix Representation Of A Feedline In The Laplace Domain 



tank. An expression of the type presented in Equation (121) is then 
assigned to each component. From this point, a minimal amount of 
matrix manipulation is then required to arrive at an expression 
similar to Equation (125), hut which reflects the particular configura- 
tion under consideration. The flow of subsequent calculations in the 
program is governed by the following inputs: 

(1) A coded set of integers describing the type of elements 
in the order in which they appear in the line, 

(2) The number of Bi's that are to be summed to arrive at 
matrix B_ and 

(3) The number and type of matrices that are contained in 
each B_i . 

Using the codes in Item (1), additional data, relevant to each compo- 
nent type, are read in and associated with a specific component lo- 
cation in the line. Categorically, this data can include component 
lengths, line radii, friction factors, bubble radii, spring constants, 
damping constants, liquid-vapor mass ratios, etc. The complete 
input package also includes information that pertains to the line in 
general, e. g. , input and terminal impedances, thermodynamic and 
fluid dynamic properties of the propellant, mean flow velocity, the 
elastic and geometric properties of the conduit wall, and the fre- 
quency band over which the response is to be determined. 
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VII. EXPERIMENTAL VERIFICATION OF 
THE FEEDLINE MODEL 


An experimental test 'program has been conducted to determine 
the validity of the analytical model and corresponding computer program 
for a liquid rocket propellant feedline. A basic test apparatus, which 
simulated in an elementary way a typical liquid rocket feed system, was 
constructed and assembled; this facility was sized so that the test 
parameters spanned the ranges of interest. From initial baseline tests, 
the system complexity was expanded by introducing various effects, and 
experimental responses were obtained for comparison with the com- 
puter analysis. 


VII. 1 Experimental Feedline Test Facility 

The basic objectives in the design of the experimental feedline 
test facility were to approximately simulate the pressure response 
(amplitude and frequency) generated by a typical feedline excited by 
various means, while using water as the internal fluid rather than an 
actual rocket propellant or LOX. One feedline perturbation technique 
is illustrated in Figure 21 in which the feedline is externally excited by 
introducing a dynamic volume flow perturbation superposed on the steady 
mean flow. This technique has been used by other researchers l > 19 > 20 
to obtain the dynamic response of specific feedline configurations, and 
it was chosen as the excitation scheme for the initial tests conducted in 
this experimental program due to its relative controlability. 

Dimensional analysis of a problem involving a typical LOX feed- 
line with a finite terminal resistance and excited by a flow pulsar reveal- 
ed that the following similarity parameters should be observed for exact 
simulation by the experimental facility: 



Pi 

(126) 

C °« Pm 

c t p 

(130) 

Qm 

" Z t Q 

Z tB Qm 

Z t Q 

Q n 

Qi 
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(127) 
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Figure 21. Illustration of Feedline Excited by a Pulsar 
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where m = model. 


Ideally, it was initially desired to simulate a typical LOX feedline 
approximately 30 ft long and 8 in. in diameter with water as the fluid 
medium in a line 3 in. in diameter. The 3 -in. diameter line was chosen 
based on available pumping equipment and power requirements. This 
allowed for the generation of a turbulent mean flow at the test conditions, 
which was necessary for verification of the turbulent flow propagation 
operator model developed for the computer program. The actual length 
of the line was 27 ft, and this dimension was chosen to obtain a first 
longitudinal mode resonant frequency with water in a rigid line of approxi- 
mately 45 Hz. This provided a low enough resonant frequency such that 
the hydraulic actuator which drove the side branch pulsar piston operated 
with a suitable amplitude so as to produce a reasonable (measurable) 
pressure response slightly past the first anti- resonance. 


From the scaling laws presented above, the response of a LOX- 
filled line of the same relative geometry as that of the model can be easily 
obtained. The primary frequency responses obtained during the test 
were the terminal impedance amplitudes, defined as the dynamic pres- 
sure perturbation amplitude, P 4 , at the line terminus divided by the 
input dynamic volume flow perturbation amplitude, Q d . 

Dividing the similarity parameter in Equation (136) by Equation 
(133)* the corresponding LrOX line pulsar impedance amplitude can be 
obtained from the experimental model amplitude as 


Qd 




(137) 


Combining Equations (129) and (130), the ratio of the terminal re- 
sistances required for the simultaneous simulation of viscous effects 
(Reynolds number) and compressibility (acoustic effects) becomes 



P 3 Si kind 3 
p B 3 c 0 2 b |ad 3 


(138) 
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while the corresponding frequencies are determined from 


oj - 




( 139 ) 


From the remaining similarity or scaling laws, all the parameters that 
would exist in the simulated LOX feedline can be computed; each vari- 
able in the model has a corresponding value in the simulated LOX feed- 
line case. 


A schematic of the feedline flow facility used for the basic tests 
is presented in Figure 22. Shown are the relative locations of the lower 
reservoir, pump, hydraulic motor, upper reservoir, supply line and 
feedline, the pulsar piston, hydraulic actuator and terminal impedance 
or resistance. The photograph in Figure 23 reveals the actual relative 
sizes of each of the lower components, including the tower and .support 
structure. The flow was produced by a five- stage, deep- well type tur- 
bine pump, capable of producing 600 gallons /minute with a total dynamic 
head of 460 ft at 3500 rpm. The pump was powered by a Vickers fixed- 
displacement, piston-type hydraulic motor rated at 100' HP at 3600 rpm. 
The side branch pulsar piston was driven by a MOOG hydraulic servo 
valve and actuator assembly with a static force capability of 10, 000 lb 
force decreasing to 6, 000 lb at 100 Hz. The maximum stroke of this 
unit was 2. 84 inches. The piston area was designed large enough to 
produce a dynamic pressure amplitude equal to 50 psi {100 psi peak-to- 
peak) with the use of a relatively low terminal resistance, Zt = 1. 9 X 10 5 
lb- sec/ft 5 . 


The transfer function relating the pressure response at the line 
terminal to the input dynamic flow perturbation, Q d , for the line con- 
figuration used in the base line tests is given by 

fit ZcA" 1 sinh (YL) • ... 

Q d ” (Z c /AZ t ) sinh (yL) + cosh (YL) 1 ; 

For an in viscid model, the terminal impedance amplitude at resonance is 
simply equal to the value of the terminal resistance. 


Qd 


= z fc 


(141) 


and from this expression, the required volume flow perturbation ampli- 
tude was found to be Q d = 0.076 ft s /sec. It was desired to conduct a 
given test with constant input perturbation amplitude over the frequency 
range; therefore, the piston stroke required was given by 
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Figure 22. Flow Facility for Current Experimental Program 
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Figure 23. Lower Portion of Feedline Facility 
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(142) 


X„ 


Qd 

2nf A p 


where f is the input or excitation frequency and A p is the piston frontal 
area. Figure 24 presents the hydraulic actuator stroke limit, which is 
based on a maximum servo valve flow amplitude of 44. 8 cubic inches per 
second (peak-to-peak) at 100 Hz and an internal actuator piston area 
equal to 3. 3 inches 2 . Also shown is the stroke required versus frequency 
for the pulsar piston with diameter equal to 4 inches. 


The upper reservoir was fabricated from a 77-gallon stainless 
steel tank, rated at 400 psi. For all tests conducted, an inlet impedance 
approximately equal to zero was provided by trapping air in the top of the 
tank so as to provide a liquid-gas interface. The base of the reservoir 
support structure was 30 ft above ground level, and was supported by a 
triangular tower composed of three sections, each 10 ft long. The actual 
length of each line was 25 ft, the remaining 2 ft being 1 ft for the pulsar 
"tee" section and the other length connected to the tank outlet. Figure 25 
presents a photograph showing a close-up view of the pulsar "tee" section, 
the hydraulic actuator, and the pulsar piston cylinder. The piston was 
sealed with an "O" ring, and lubrication was provided from the rear by 
filling the remainder of the cylinder bore with hydraulic oil, which was 
displaced periodically into the attached accumulator. 


Instrumentation 


A simplified schematic of the basic instrumentation required for 
the experimental testing is shown in Figure 26. The flow of hydraulic 
pressure into the actuator was controlled by the servo controller. The 
servo control system consisted basically of a DC "set point" input which 
set the hydraulic actuator at the mid-range of its stroke, plus a sinu- 
soidal input to generate the flow perturbations. The feedback from the 
piston was obtained by means of a velocity transducer. The servo con- 
troller electronic circuit is shown in Figure 27. The servo controller/ 
valve was driven by an electronic oscillator, from which the input fre- 
quency and amplitude could be manually controlled. The input flow 
perturbation amplitude was measured with a Sanborn linear velocity 
transducer, which had a sensitivity of 94 millivolts /inch/ second and a 
maximum stroke of 1. 0 inch. From the known piston area, the flow 
perturbation amplitude, Qd was calculated in cubic ft/second. The 
flow perturbation amplitude was displayed on the oscilloscope and con- 
trolled manually at each frequency setting since the control loop was 
open. The velocity transducer can be seen in Figure 25. 

The dynamic pressure response was measured through the em- 
ployment of Kistler 603A piezoelectric quartz dynamic pressure 
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ACTUATOR STROKE (INCHES) 



Figure 24 . Stroke Required to Produce a Flow Perturbation Amplitude 
of .076 fr/sec 
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Figure 25. Hydraulic Actuator and Pulsar Assembly 
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Figure 27. Servo Contrail 


transducers; the charge signals were converted to voltage readings with 
Kistler Model 504D electrostatic charge amplifiers. Special ultra-low 
noise cables were used to interconnect the transducers and amplifiers. 
These transducers were small with . 22 mm diameter faces, and were 
mounted flush with the internal feedline wall. These transducers had a 
resolution of 0. 01 psi and were acceleration- compensated to null out 
spurious vibrations. Each transducer was calibrated with the charge 
amplifier to determine its actual sensitivity in psi/volt. 

It was discovered during early testing that the dynamic flow per- 
turbation inputs were somewhat distorted from pure sinusoidal wave- 
forms; this distortion was produced by mechanical friction along the 
hydraulic actuator shaft and pulsar piston. These input distortions not 
only made reading the correct flow perturbation amplitude difficult, but 
also produced distortions in the pressure signals and made measurements 
of its phase lag doubtful. A two- channel SKL 308A variable electronic 
filter was used to clarify both the input and pressure response signals. 

The low-pass cutoff frequency for each channel was set at 1. 5 times the 
excitation or input frequency. In this manner, the phase lags created by 
the filters were equal, keeping the phase lag between the input signal and 
the pressure response equal to its unfiltered value. 

The dynamic pressure and flow perturbation signals were then 
recorded on a CEC Type 5-116-P4-14 oscillograph recorder, which was 
operated at various paper speeds ranging from 8 to 64 inches/ second, 
depending upon the input frequency. Each channel was calibrated so as 
to obtain the inches deflection/ volt input, and from the known pressure 
and velocity transducer sensitivities, the pressure and velocity sensiti- 
vities were obtained as psi/ inch and cubic ft per second/inch, respec- 
tively. Figure 28 presents typical raw oscillograph data obtained, with 
and without the electronic filters in the system. Also shown is the method 
used to determine the phase angle, 0, from the raw data records. 

VII.2 Base Line Test Results 


The feedline configuration used for the base line tests is presented 
in Figure 22. The initial tests were conducted with the termination of 
the line blocked. Also, the initial feedline employed was fabricated from 
thick-walled (Schedule 40S) type 304W stainless steel pipe, 3 in. nominal 
pipe size. With this material and geometry, the effects of the radial wall 
compliance were negligible or slight. By blocking the line terminal, 
possible effects of a mean steady flow component were eliminated from 
consideration, and the problem reduced to that of a closed hydraulic con- 
duit connected to a large reservoir at one end and excited by a flow per- 
turbation at the closed end. This simplified problem allowed for an 
evaluation of the instrumentation and testing techniques applied since 
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for this case, excellent agreement between theory and experiment should 
exist, as was the result obtained by Gerlach in Reference 18. This was 
the most elementary configuration tested. For this test, the pressure 
amplitude and its phase angle with the dynamic flow perturbation input 
was recorded at the line terminal and also in the reservoir. The pres- 
sure amplitudes generated in the upper reservoir, even at line resonance, 
were insignificant, indicating that the assumption of zero impedance at 
the line inlet was reasonably correct. 

Figure 29 presents the frequency response obtained for a stain- 
less steel line 27 ft long and 3. 068 in. in diameter, with a wall thickness 
equal to 0. 216 in. , and excited by the side branch flow pulsar. The 
pulsar flow perturbation double amplitude was set at Qa = 0. 007 ft 3 / sec, 
while the excitation frequency was increased from 5 to 100 Hz. Since the 
flow perturbation amplitude proved difficult to control at this exact ampli- 
tude for all frequencies, the resultant pressure amplitude was divided by 
the input flow perturbation amplitude to produce a dynamic terminal im- 
pedance amplitude, P T /Q d . The linearity of the system thus allows the 
pressure response at any frequency and input amplitude to be obtained 
from Figure 29. The terminal pressure response phase lag was also 
obtained and is plotted as a function of input frequency in Figure 30. The 
data were obtained by pressurizing the liquid to 100 psig by pumping the 
water into a "dead" head; 100 psig was chosen to eliminate any possible 
cavitation problems which might have been developed by the large ampli- 
tude dynamic pressures at resonance. Figure 29 indicates excellent 
agreement between theory and experiment, except there appears to be 
more damping at resonance than predicted by the analytical model with 
an infinite terminal resistance. The apparent disagreement in phase 
angle near resonance can be attributed to the poor measurements of the 
phase relationships from the data traces, as these data were recorded 
before it was realized that electronic filtering of the signals was re- 
quired as discussed earlier in Section VII. 1. 

Effect of the Terminal Resistance 


To determine the effect of the terminal resistance, several tests 
were conducted with various values of the terminal impedance, AF/Q. 
Tests were accomplished with both linear and nonlinear resistance ele- 
ments; the resistance for small perturbations from the steady mean 
value was as sinned linear in the analytical model, and the validity of this 
assumption was to be verified. This is important for a liquid propellant 
feed system since the pump and injectors being simulated by the resis- 
tance element actually have nonlinear steady flow characteristics which 
are generally linearized for modeling purposes. 
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Figure 29 . Frequency Response of a Blocked Stainless Steel Line 
Excited by a Side Branch Pulsar 
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Figure 30 . Phase Lag for a Blocked Stainless Steel Line Excited by a Side 
Branch Pulsar 
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By using various degrees of resistance in the experimental tests, 
various mean flow velocities were also established. In most cases these 
flow velocities were large enough that turbulent flow was established in 
the system, which was required to determine the, success of the modified . 
propagation operator for a turbulent mean flow discussed in Section IL 3. 


Tests with Linear Flow Resistance 


A considerable effort was directed toward the development of a 
suitable linear flow resistance element for the line terminus. The basic 
requirement for a linear relationship between pressure drop and flow 
rate through a hydraulic element was that the flow be laminar in the local 
passages of that particular component. Laminar flow required a local 
Reynolds number less than approximately 2000 for internal flows. From 
the definition of Reynolds number. 



(143) 


a low Reynolds number was obtained for a given fluid and velocity by de- 
creasing the passage diameter, d. 


Initial attempts at constructing a linear- resistance was through- the 
use of fine mesh wire screen (325 mesh). This proved successful for low 
values of resistance, AP/Q. However, as the flow velocity and pressure 
drop increased, the relation between pressure and flow approached that 
of an ordinary orifice, AP = kQ 2 . A suitable linear resistance was finally 
obtained by packing assorted diameter glass beads between two 325 mesh 
wire screens and into a cylinder 3 inches in diameter and 4 inches long. 
This length was selected based on the experimental data collected in 
Figure 31 for several attempts at developing a true linear pressure-flow 
relationship. Figure 32 presents the actual pressure drop versus volume 
flow measured for this configuration, in which the bead diameters ranged 
from 0. 0117 in, to 0. 0029 in. , providing very small diameter flow 
passages. The 4-in. thickness provided a much more linear flow-pressure 
behavior than any of the other configurations shown, with the pressure drop 
and flow being related by AP = kQ^*^. The trend in the summary calibra- 
tion of the various "linear" resistance elements tested indicates a true 
linear resistor might be approached or obtained but would require a con- 
siderably long or thick resistor. Increasing the length also increased the 
magnitude of the resistance, requiring an extremely high pressure drop 
for a given flow rate. For these reasons, the 4-in. thick resistance 
element was selected. 


The tests conducted with this linear resistance element were quite 
similar to those described for the base line test with an infinite terminal 
resistance. Figure 33 presents the frequency response obtained for an 
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Figure 32. Calibration of "Linear" Resistor 
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FREQUENCY, HERTZ 


Figure 33. Frequency Response of an Aluminum Line with a Linear 
Terminal Impedance, Excited By a Side Branch Pulsar 
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aluminum, line 27 ft long and with a diameter of 3 in. , excited by the side 
branch pulsar located just above the line terminal. The wall thickness 
was 0. 050 in. and the terminal impedance, Z t , was 9. 75 X 10 s lb-sec/ft B 
producing a mean flow rate Q = 0. 0108 ft 3 /sec at 80 psid pressure drop 
through the resistor. The corresponding phase lag versus frequency is 
shown in Figure 34. The predicted amplitude, ‘phase, and resonant fre- 
quencies agree quite well with the experimental data for the first and 
second longitudinal modes. 

Tests with Flow - Nonlinear Resistance 


Tests were conducted with the feedline configuration described 
above for nonlinear type resistors at the line terminal. These resistors 
were simple orifice plates; one designed to produce a terminal impe- 
dance value, Z t , near that of the already described linear resistor and 
another had a nominal impedance much lower. The resultant terminal 
pressure responses versus excitation frequency were obtained with 
7/16-in. and 1/4-in. diameter orifices installed in the 3-in. internal 
diameter feedline. The discharge coefficient, k d , for the resultant 
diameter ratio, 8, and Reynolds number was obtained from standard 
tables^* and checked by a direct calibration while installed at its posi- 
tion in the flow facility. With a knowledge of the discharge coefficient, 
the pressure drop versus mean flow rate was calculated from 

Q = k d A.^' (144) 

where A 0 was the cross-sectional area of the orifice. From the defini- 
tion of the terminal impedance, Ap/ Q, the value of the resistance at 
the operating point was obtained by writing Equation (144) as AP = f(Q 3 ) 
and differentiating to obtain the slope of the pressure drop versus flow 
curve. Since the orifice had a nonlinear pressure-flow characteristic, 
the magnitude of the resistance varied with the mean flow rate. The 
nominal value of the resistance or terminal impedance was therefore 
computed from 


7 d(AP) pQ 
* " dQ " (k d A 0 ) 2 


(145) 


where Q was the mean flow rate at the operating or test conditions. 

The frequency response of the 2 7 -ft long aluminum feedline des- 
cribed earlier was obtained for two nonlinear resistance elements 
(orifices). The 7/ 16 -in. orifice provided a nominal terminal impedance 
of Z T = 3.91 X 10 s lb-sec/ft 5 at a pressure drop of 100 psid; the ter- 
minal pressure response is presented in Figure 35 and the phase lag be- 
tween the pressure and input flow perturbations is shown in Figure 36. 
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Figure 34. Pressure Phase Lag of an Aluminum Line with a Linear 
Terminal Impedance- 
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Figure 36. Phase Angle vs Frequency for a Thin Walied Aluminum Feed- 
line Terminated by an Orifice 
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For all data points, tlie amplitude of the dynamic flow perturbations 
created by the side branch pulsar was maintained at approximately 
Qd = 0. 024 ft 3 / sec. The second orifice had a diameter of 1/4-in. , pro- 
ducing a much higher terminal resistance at any given mean pressure 
drop measured across the orifice. Tests were conducted on the same 
feedline as used for the 7/16-in. orifice plate with a steady state 
operating pressure drop of 80 psid with the 1/4-in. diameter orifice. 

Figures 35 through 38 indicate fairly good agreement with the 
analytical computer model response, plotted as the solid curve in all 
figures. At the line longitudinal resonant frequencies there appears to 
be greater damping in the experimental tests than predicted by the 
computer simulation of the same problems. Later tests conducted with 
an aluminum line with a smaller wall thickness indicate that this dis- 
crepancy cannot be necessarily attributed to the nonlinearity of the 
terminal resistance. 

An aluminum line with a 3 -in. internal diameter, 27 ft long, anc 
a wall thickness, h = 0. 022 in. , was tested with the 7/16-in. diameter 
orifice employed as the terminal impedance element. The steady flow 
pressure drop at the line terminus was 80 psid, producing an operating 
or nominal terminal impedance, Zt = 3.47 X 10 5 lb-sec/ft 6 . The mea- 
sured frequency response of the terminal pres sure /input flow perturba- 
tion is plotted in Figure 39 along with the analytical model prediction. 

The experimentally measured phase lag is also given in Figure 40 along 
with the corresponding theoretical response. It should be emphasized 
that here the nonlinear resistance was again linearized for the computer 
computation and excellent agreement with theory resulted. The lower 
amplitudes of the response at resonance in Figures 35 and 37 are pro- 
bably the results of experimental inaccuracies, as it was noted during 
several tests that the pressure response amplitudes. near resonance were 
difficult to repeat consecutively. Later tests with nonlinear resistance 
elements also produced very good agreement with theory for all cases 
where the dynamic flow perturbation amplitudes were considerably less 
than the mean flow rate. 

VII. 3 Di s tr ibute d C omplian c e s 

Two types of distributed compliances were discussed in Section III, 
and analytical models were developed for their influence on the response 
of liquid propellant feedlines. Both effects have been examined experi- 
mentally; First, distributed wall compliance effects were examined by 
using three different feedlines of the same length and internal diameter 
but with a different wall thickness, h, and in one case constructed of a 
different material. Second, the effects of entrained gases were experi- 
mentally examined by keeping the geometric and elastic parameters of the 
feedline constant and varying the gas mass ratio, $. 
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Figure 37. Frequency Response of a Thin Wailed Aluminum 
Feedline Terminated by an Orifice 





PHASE ANGLE, 



Figure 38. Phase Angle for a Thin Walled Aluminum Feedline Terminated 
by an Orifice 
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Figure 39 . Frequency Response of a Thin Walled Aluminum Line 
Terminated by an Orifice 
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Figure 40 . Phase Angle for a Thin Walled Aluminum Feedline Terminated 
by an Orifice 
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Distributed Wall Compliances 


In general, the flexibility of the feedline wall will produce two 
possible effects which may be of concern. First, the axial wall stiff- 
ness effect would permit real wave propagation of higher order modes 
at low frequency by virtue of a coupling between the fluid and axial wall 
modes. However, it was found in Section HI. 3 that for all practical 
feedline applications this effect of the axial wall stiffness on the wave 
propagation can be neglected, since the amount of energy which is fed 
into the higher order modes at the line termination is small. Second, 
the wall compliance reduces the phase velocity, and increases the spatial 
attenuation for the longitudinal fluid wave propagation mode. The classic 
Korteweg correction to the phase velocity was incorporated into the ana- 
lytical model to account for the wall compliance; the correction being 


c 


Co 



Po si d \ 2 
Et h j 


(146) 


The validity of this correction has been determined from its effect on the 
resonant frequency of the feedline, since the first longitudinal mode fre- 
quency is 



and a shift in resonant frequency is a good measure of the phase velocity 
change, providing the lengths are equal. 

For feedlines of equal lengths, L, and diameters, d, and with 
the same internal fluid, the phase velocity and resonant frequency become 
a function of the elastic parameter, E t h. 

The responses of three feedlines terminated by an infinite resis- 
tance (blocked line) were obtained for the 304 stainless steel line with a 
wall thickness, h = 0. 216 in. and two aluminum lines of the same length 
and internal diamter, but with wall thicknesses, h = 0. 050 and 0. 022 inch. 
The response of the relatively stiff stainless line was presented in 
Figure 29, while the results for the two aluminum lines are plotted in 
Figures 41 through 44. In all cases, the excitation was provided by the 
side branch pulsar piston. A summary plot. Figure 45, reveals the 
accuracy and validity of the model for distributed wall compliances. 

Shown on this plot are the measured and calculated values of the first 
mode resonant frequencies for the three feedlines tested. The resonant 
frequencies are plotted as a function of the elastic parameter, E t h. The 
results indicate that the Korteweg correction to the phase velocity is a 
valid means for modeling the effects of the wall radial compliance for 
thin- walled feedlines. 
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Figure 41. Frequency Response of Thin-Walied Aluminum Feedline 
With a Blocked Line 


90 



PHASE ANGLE, <D , (DEGREES) 



INPUT FREQUENCY, HERTZ 


Figure 42. Phase Angle vs Frequency for a Blocked, Thin-Walled, 
Aluminum Feedline 
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Figure 44 . Phase Angle for a Thin-Walled Blocked Aluminum Line 
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Tests were also conducted to verify the use of the Korteweg 
correction to the phase velocity for practical applications, where the 
terminal resistance or impedance has a finite value and where any in- 
creased damping due to wall elasticity (which has been neglected in the 
analytical model) might show up. The response of an aluminum feed- 
line with a wall thickness, h = 0. 050 in. , and terminated with a non- 
linear terminal impedance element operating such that the nominal 
terminal impedance value, Zt = 3. 91 X 10 s lb-sec/ft 5 , is presented 
in Figure 46. The corresponding phase relation between the terminal 
pressure perturbations and the input dynamic flow perturbations is 
plotted in Figure 47. Also shown in Figure 46 is the computed theo- 
retical response of the same configuration except that the feedline was 
assumed to be infinitely rigid; it is obvious that a correction is indeed 
required to properly model wall elastic effects. The results of Figure 46, 
where the terminal impedance was a finite value, justify the conclusion 
that the increased spatial attenuation for a line with an elastic wall is 
negligible over the frequency range of most practical feedline applica- 
tions. 


Entrained Gases 


The net effect of entrained gases in the propellant or simulation 
fluid has been found to decrease the phase velocity and also increase the 
spatial attenuation by introducing added damping. In determining the 
effect of a uniform distribution of gas bubbles on the frequency response 
of feedlines, it was found that a very accurate means of measuring the 
mass of gas present in the system at the time of testing was necessary. 
Nitrogen gas was chosen for these tests, and its theoretical effect on the 
reduction of the speed of sound in the water-gas mixture was computed 
for several mixture pressures and is shown in Figure 48. From these 
computations it was evident that only a small mass of entrained gas would 
drastically vary the frequency response of a given feedline requiring a 
sensitive method for the measurement of the mass ratio. An injector was 
also designed and fabricated to generate a uniform distribution of small 
nitrogen gas bubbles assumed by the analytical model. 

The mass ratio, of the liquid-gas mixture was determined 

from measurements made of nitrogen gas mass flow rate into the feedline 
system. The gas bubbles were injected at the upper end of the feedline, 
just below the tank or upper reservoir. The mixture of nitrogen gas 
bubbles and water were then swept downstream with the mean flow at the 
rate 


v B = v 0 - v BR (148) 

where v Q was the velocity of the water in the line and vsr was the bubble 
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Figure 46. Frequency Response of an Aluminum Line Excited by a Side 
Branch Pulsar 







Figure 47. Phase Angle vs Frequency for a Line Excited by a Side 
Branch Pulsar 
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rise velocity. In all cases, the bubble rate of rise had to be less than the 
mean water velocity; the gas flow rate measurements would have been 
meaningless if an unknown loss into the upper reservoir occurred. By 
making the bubbles small, their rates of rise were well below that of the 
water velocity at the test conditions. A Cole-Parmer low volume flow 
meter which had a range of 10 to 1900 cubic centimeters per minute was 
used to measure the volume flow rate; the pressure at the meter was 
recorded and used to calculate the gas density at the meter, from which 
the nitrogen gas mass flow rate into the feedline system was obtained. 

With the feedline facility being operated "open loop", the total 
mass of flow into the system at the injector was assumed to pass down- 
ward through the feedline; i. e. , none was assumed to escape upward into 
the reservoir. To verify the assumption regarding the bubble rise velo- 
cities, a special bubble calibration chamber was designed to measure both 
the bubble sizes being generated by the injection device and their rates of 
rise. The photograph (Figure 49) shows the calibration tube, constructed 
of acrylic tubing (plexiglass) and partially filled with water. The chamber 
was pressurized to value approximating those at the feedline entrance 
(tank exit) which would occur during tests. For calibration purposes, the 
gas injector was installed at the bottom of the tube; and the nitrogen gas 
was injected at approximately the rates required for the test mass ratios, $. 
The pressure was maintained constant by allowing the gas to bleed from the 
upper volume at the same rate as it was being injected. Figure 50 shows 
a photograph taken close-up, with a scale mounted on the chamber wall 
from which measurements were made of the bubble sizes. From these 
visual observations, it was determined that the bubble distribution was 
nearly homogeneous. The bubble rates of rise were obtained by dividing 
a known distance marked on the calibration chamber wall by the time re- 
quired for a selected group of bubbles to traverse that distance. The 
bubble rise velocities varied from 0. 2 to 0. 4 ft/ sec, depending on the 
bubble size. 

A close-up photograph of the bubble injector is presented in 
Figure 51. This injector was designed to bolt between the flange at the 
tank exit and the upper feedline flange. The nitrogen gas entered the 
injector at two ports located 180° apart, and passed around a channel cut 
circumferentially around the injector flange. The bubbles were finally 
created by passing the gas through small diameter tubes, fabricated from 
hypodermic tubing and drawn down by heating and stretching one end until 
the diameter was of the order of 0. 001 inch. Other researchers, 

Carstensen and Foldy^, and Fox, Curley, and Larson^, performed 
tests requiring a uniform distribution of gas bubbles in a liquid; it was 
from their experience that ideas were drawn for the development of this 
gas injector. 
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Figure 49. Calibration Tube for Bubble Measurements 
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Figure 50. Close-Up of Bubble Sizes and Distribution 
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Figure 51. Nitrogen Gas Bubble Injector 
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Figure 53. Phase Angle vs Frequency for an Aluminum Line with a 
Water-Nitrogen Gas Mixture 
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Figure 54. Frequency Response of an Aluminum Feedflne with 
a Water-Nitrogen Gas Mixture 
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Figure 55. Phase Angle vs Frequency for an Aluminum Line with a 
Water-Nitrogen Gas Mixture 



than for the § = 4. 31 X 10~ 8 and $ = 0 cases, indicating a possible in- 
creased damping with mass ratio, $. It should be noted, however, that 
the difference in amplitude between theory and experiment at resonance 
is really no greater than that observed without bubbles, as evidenced by 
Figures- 35 and- 37,' The' difference could' be only' a result of experi- 
mental inaccuracies in measuring the correct amplitude at resonance. 

From Figure 48 it can be observed that near $ = 10~ 5 , the rate 
of change of the mixture speed of sound with mass ratio is near its maxi- 
mum; hence, a small error in the measurement of mass ratio, $, could 
have produced the apparent discrepancy between the theoretical and ex- 
perimental resonant frequencies for this configuration and shown in 
Figure 54. The experimental data indicate that the mass ratio was 
larger than that measured and thought present. Data were also recorded 
for a measured mass ratio, § = 5.61 X 10 -6 , which was the greatest 
mass ratio produced for any of the experiments conducted. The results 
are shown in Figures 56 and 57. Again, the data showed a lower reso- 
nant frequency for both the first and second longitudinal modes than the 
theoretical predictions, although in this case the resonant frequency was 
off by only 1. 5 Hz for the first mode and 5 Hz for the second peak. The 
pressure/flow perturbation amplitudes at resonance were considerably 
damped with this much greater volume of distributed nitrogen gas in the 
system. The theoretical response also indicated more dissipation of 
energy at resonance, but the degree of damping was significantly less 
than that which was obviously present. It is possible that further re- 
finements should be made to the analysis to reflect this observed in- 
crease in attenuation for the larger mass ratios of gas to liquid. This 
increased attenuation with increasing gas/liquid mass ratio, i, suggests 
that injection of inert gases, as suggested in earlier research^, might 
be a very suitable method of POGO suppression. Measurements made and 
reported in Reference 24 indicate that for reasonable volume fractions 
(volume gas/volume of liquid) the thrust produced is hardly affected. 

VH.4 Local Compliances 

The next phase of the experimental test program conducted for 
verification of the analytical model and computer program involved gen- 
erating local compliances or compliance -type elements and experimentally 
measuring their effects on the frequency response of feedlines. More 
specifically, the. validity of the model used to simulate "complex" side 
branches and localized gas or vapor bubbles was established. The term 
complex side branch refers to any type branch in which there exists a 
fluid resistance and/or inertance as well as a compliant-type element. 

The analytical models for both categories of local compliances were 
discussed in Section IV of this report. 
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Figure 56. Frequency Response of An Aluminum Feedline 

with A Water - N 2 Mixture 
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Complex Side Branch 


The types of side elements typical of rocket feet systems (accumu-' 
lators and partially gas-filled pressure-sensing lines) behave as local 
compliance -like elements. The experimental arrangement used to verify 
the computer simulation (analytical model) is shown in Figure 58. The 
side branch was located 4 feet above the line terminal; the length of the 
branch was 10. 2 feet and the diameter was 1.4 inches, with a wall thick- 
ness equal to 0. 049 inch. A short (4-inch long) plexiglass tube was 
attached vertically at the end of the branch; this allowed for visual obser- 
vation and measurement of the volume of trapped gas at the branch termi- 
nus. Tests were performed to obtain response curves (amplitude ratio 
and phase angle) for the main feedline, with the system excited by the 
side branch pulsar. In addition, dynamic pressures were recorded at 
the end of the side branch for several of the tests. 

The feedline geometry described above was tested with water as 
the fluid media, and various volumes of gas (air) were trapped at the side 
branch terminus. Figure 59 presents the frequency response obtained 
with a terminal impedance, Zt = 3. 91 X 10 s lb-sec/ft g , and with a branch 
terminal volume equal to 0. 769 in. 3 . It appears that the measured reso- 
nant frequencies have been shifted (lowered) by the amount predicted by 
theory, but the damping effect created by the side branch is somewhat 
greater than that predicted by the analytical model. The phase relation 
between the terminal pressure perturbations and the input flow perturba- 
tions are plotted in Figure 60. 


Figures 61 and 62 present the measured pressure-to-flow pertur- 
bation ratio amplitudes and their corresponding phase relation for the 
same line geometry with a larger air bubble volume (V 0 = 2.67 in. 3 ) at 
the side branch terminus. These data also generally agree with the 
theoretical results, also plotted in the same figures as solid lines, ex- 
cept the experimentally determined amplitudes at resonance are some- 
what less than expected. The- effects of a further increase in bubble size 
to a bubble volume equal to 5. 25 in. 3 with the same terminal impedance 
and feedline geometry are presented in Figures 63 and 64, In general, 
these data represent the poorest agreement with theory of all the data 
collected for a feedline with a complex side branch. 

Close examination of Figures 59 through 64 reveal that the net 
effect for an accumulator-type complex side branch element, terminated 
by a trapped volume of gas, is to shift (lower) the resonant frequency in 
proportion to the trapped gas volume and to also produce increased 
damping due to the resistance in the branch line. Perhaps an improved 
model of side branch elements would be to treat the branch line as a 
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Figure 58. Geometry of Feedline Excited by a Pulsar with a Complex 
Side Branch 
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Figure 60 . Phase Angle vs Frequency for a Line with a Side Branch 
Element with Trapped Air 
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Figure 62 . Phase Angle vs Frequency for a Line with a Side Branch with 
Trapped Air 
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Figure 63. Frequency Response of a Feedline with a Side Branch 
with Trapped Air 






Figure 64. Phase Angle vs Frequency for a Line with a Side Branch 
Containing Trapped Air 
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distributed parameter system rather than as a lumped element as des- 
cribed in Section IV. 

The above described tests experimentally measured the effects of 
an accumulator-type side branch on the frequency response of a simu- 
lated feedline. A separate test was also conducted to simulate a pressure- 
sensing line partially filled with gas. This side branch element was of a 
much smaller diameter than that used in the previous tests; the diameter 
was 1/4 inch, and the overall branch length was 89 inches. This branch 
was connected to the primary feedline at the line terminus, 1 inch above 
the resistance element and at the same longitudinal station as the pres- 
sure transducer used to monitor the terminal pressure perturbation ampli- 
tude. The branch line was curved slightly upward and was partially filled 
with air trapped between the liquid column and the branch terminal. A 
pressure transducer was also mounted at the end of the branch line to 
measure the dynamic pressure at the branch terminal in order that it 
might be compared with the actual pressure in the feedline at the same 
station. The volume of trapped air was measured to be 1. 96 in. 3 . The 
perturbation pressure/flow amplitude and phase versus excitation fre- 
quency are shown in Figures 65 and 66, respectively, for this configura- 
tion. The sharp, low-frequency peak at 4. 6 Hz predicted by theory was 
not detected or measured with the original test conducted^ since no data 
were collected below 5 Hz. After the computer simulation of the same 
problem was obtained, the response at 4. 6 Hz was observed on the com- 
puter print-out, and the test was set up and re-conducted and data were 
obtained from 1 Hz to 7 Hz. Excluding the low-frequency response, the 
measured data shown in Figures 65 and 66 agree extremely well with the 
analytical model. 

The experimentally determined transfer function relating the pres- 
sure amplitude measured at the branch terminus to that actually existing 
in the primary line at the same longitudinal station (terminal) is plotted 
in Figure 67. This problem simulates a remotely located pressure sensor, 
connected by a partially gas-filled pressure-sensing transmission line; 
it is obvious that care should always be taken in determining the transfer 
function relating the pressure amplitude at the remote sensor to that at 
the actual feedline in order that the correct pressure amplitudes in the 
feedline are determined. 

Large Localized Gas Volumes (Bubbles) 

Section IV. 1 presented the analytical model developed for localized 
compliances, such as a large pocket of trapped gas or cavitation bubbles. 
For this particular model, it was assumed that a cavitation region could 
be approximated by a large spherical bubble, and that the condensation and 
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Figure 65. Frequency Response of a Feedline with a Side Branch Containing 
Trapped Air and Excited by a Flow Pulsar 





Figure 66. Phase Angle for a Feedline with a Side Branch Containing 
Trapped Air and Excited by a Flow Pulsar 
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Figure 67. Experimental Transfer Function Relating Side Branch Pressure 
Amplitude to the Feedline Terminal Pressure 
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evaporization rates were such that a single -component, two-phase mixture 
of a liquid and its vapor could be treated as a constant quality process; 
that is, there is no change of state during successive expansions or con- 
tractions of the cavitation region. For a two-component mixture, such as 
BOX and a large localized region or bubble of captured helium ullage gas, 
the above 'assumption is -unquestionably valid. 

Tests were conducted with the experimental feedline test facility to 
measure the effect of a local compliance located at an arbitrary position 
longitudinally along the feedline. Since the model used in the computer 
simulation assumed a spherical gas bubble of constant mass (no phase 
changes), it was desired to create a spherical bubble and keep it stationary 
in the feedline at any designated longitudinal station. This was accom- 
plished by using a nearly spherical rubber balloon; a schematic showing 
how the balloon was held in place and the apparatus used to measure its 
compliance or pressure and volume is given by Figure 68. A plexiglass 
section 4 inches long was inserted into the feedline at a point just above 
the side branch pulsar, or 1 foot from the line terminal. 

In order to calculate the gas volume compliance, accurate mea- 
surements of the balloon volume and pressure were required at the test 
conditions. To accomplish these tasks, the bubble was initially filled to 
an arbitrary volume by pressurizing it with the feedline operating pres- 
sure near 1 atmosphere. The feedline was then pressurized to its desired 
internal test pressure. Valve B was then slightly opened, and the external 
pressure on the rubber balloon from the fluid forced the balloon to collapse 
until its volume was either zero or negligible. Valve B was then closed. 
The reservoir volume and the volume of the line connecting valves A and 
B, as well as the volume of the line to the pressure gage, were measured 
previously and recorded as V c . Valve A, which connected this apparatus 
to a high-pressure supply of nitrogen gas, was then opened and the known 
volume region was pressurized to approximately 10 psi over the pressure 
existing in the feedline; this pressure, P^ , was then recorded after 
valve A was closed. Valve B was then again opened and the entire system 
was allowed to come to equilibrium at pressure P 2 . From Boyles' law, 

Pj V c = P 2 V 3 (150) 

where 

v 3 = V 0 + V 0 , (151) 

and the balloon volume was calculated from 
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Figure 68. Apparatus for Generation and Measurement 

of a Local Compliance 
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This technique allowed for an accurate measurement of the local 
compliance and also provided a method for keeping the bubble at a station- 
ary longitudinal position along the feedline* Due to the surface tension in 
the balloon skin, the bubble spring rate (stiffness) could not be given simply 
by 


k = YPo/V 0 (153) 

but needed to be corrected by the addition of a wall tension correction 
factor, k * The total balloon or bubble stiffness, k, became 

k = YPo/V 0 + k' (154) 

where k was determined experimentally by measuring -AP/AV for the 
balloon under static conditions. From the measurements obtained pre- 
viously for P Q and Y 0 , the added stiffness, k # , was then determined. 
This additional stiffness was determined to be relatively small when com- 
pared with the uncorrected value, k = YP o /V 0 , at the operating pressures 
and bubble volumes used in the experimental tests. 

Tests were conducted with various values of the local compliance 
(bubble volume). An aluminum feedline 27 feet long, 3 inches in dia- 
meter, and with a wall thickness equal to 0. 050 inch was used in the 
examination of local compliance effects. The frequency response of this 
same feedline without a localized compliance was presented in Figure 46. 
Figures 69 and 70 present the pulsar impedance frequency response and 
phase angle obtained for the first test case in which the bubble volume was 
1.645 cubic inches and the local pressure, P c , was 111.5 psia. Again, 
the pulsar impedance amplitude refers to the terminal pressure perturba- 
tion divided by the input flow perturbation. The terminal resistance, Z t , 
was equal to 3.91 X 10 s lb-sec/ft s . The agreement between the analytical 
model and the experimental data is excellent; the measured amplitudes at 
resonance are only slightly below that predicted by theory, and the first 
and second mode resonant frequencies were predicted by the computer 
simulation within 1 cycle per second or less. 

Figures 71 and 72 present data obtained for a bubble whose volume 
was measured to be 4. 88 cubic inches; the pressure was measured to be 
109 psia. The bubble was located longitudinally at the same station as that 
previously tested, 1 foot above the line terminal. The terminal impedance 
was 3. 91 X 10 5 lb- sec/ft 5 ; the approximate average flow excitation ampli- 
tude was 0. 02 ft /sec. As predicted by the analytical model, increasing 
the local compliance value (bubble volume) lowered the first mode resonant 
frequency and sharpened the second mode peak. Figures 73 and 74 show 
the measured effect of still further increases in the bubble size. This 
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Figure 70. Phase Angle Between the Terminal Pressure and Pulsar Flow 
Perturbations fora Line with a Local Compliance 
(Bubble Volume = 1.645 in. 3 )* 
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Figure 72. Phase Angle Between the Terminal Pressure and Pulsar 
Flow Perturbations for a Line with a Local Compliance 
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Figu re 73. Frequency Response of a Feedline with a Local Compliance 
at the Line Terminal and Excited by a Flow Pulsar 
( Bubble Volume = 5.70 in^ ) 
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Figure 74. Phase Angle vs Frequency for a Feediine with a Loca ! 
Compliance 


o°° 




bubble was measured to be 5. 70 cubic inches at an internal pres-sure equal 
to 111.2 psia. The bubble was also located 1 foot above the terminal re- 
sistance element, and the feedline was excited by the side branch flow 
.pulsar... The amplitude at. the. first longitudinal mode. r.es.onant. frequency — 
agrees very well with the analytical solution, while the second mode reso- 
nant frequency agrees within 3.5%, but appears to be lower than the pre- 
dicted amplitude. 

Figures 75 and 76 present the pulsar impedance frequency response 
and phase angle obtained for a bubble located at a position other than at the 
end of the line. For this experiment, the localized volume was 5. 75 cubic 
inches and at an internal pressure equal to 111.5 psia. The nitrogen-gas- 
filled balloon was located 4. 55 feet above the line terminus; the geometry 
'of the feedline was exactly that used for the other three tests. Note that 
the volume, and hence compliance, of this bubble is almost exactly equal 
to that described in Figure 73. As evidenced by a close comparison of 
that data with Figure 75, the major effect of bubble longitudinal location 
is to increase the first mode resonant frequency and increase the response 
roll-off after resonance as the bubble is moved up the line toward the tank. 
The second mode response has a more gradual rate of increase as reso- 
nance is approached than did the response for the same bubble volume lo- 
cated nearer the line terminal,. 

It should also be mentioned that the curves labeled "theory' 1 in all 
the data collected for localized compliances did not include the previously 
discussed stiffness correction factor, k'; yet there is good agreement 
with' the experimental data since k' was always small compared with 
YP 0 /V„. The presence of the balloon wall between the gas and liquid also 
undoubtedly affected the natural damping of this compliant element, but 
since the bubble resonant frequencies were considerably higher than any 
of the test frequencies, the bubble dynamics were negligible, indicating 
that a more simple model can be used for a local compliance when the 
expected excitation frequencies are far below the bubble resonant fre- 
quencies . 

VII. 5 Structural- Hydraulic Coupling 

Experimental tests were conducted to determine the effects of 
structural -hydraulic coupling produced when the feedline was excited ex- 
ternally. A feedline was fabricated from aluminum pipe; this feedline had 
two 90° bends (short radii elbows), and its geometry was as is shown in 
Figure 77. 
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Figure 75. Frequency Response of a Feedline with a Local Com 
pliance Located Above the Line Terminal and 
Excited by a Flow Pulsar 
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Figure 76 . Phase Angle Between the Terminal Pressure and Pulsar 
Flow Perturbations for a Line with a Local Compliance 
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Line with Bends 


To determine the effect of bends or elbows on the frequency re- 
sponse, the feedline was installed as shown in Figure 77. For this experi- 
ment, iFwas desired^ to measure the effect~bFthe bends without external 
structural acceleration; hence, the excitation was provided by the side 
branch pulsar, driven by an hydraulic actuator. The resulting frequency 
response (pulsar impedance amplitude, P T /Q d ) and the phase relationship 
between the input flow perturbations and the terminal pressure perturba- 
tions are presented in Figures 78 and 79. The terminal resistance was 
created by a 7 / 16 -inch- diameter orifice which produced, when operated 
at a steady state flow pressure drop of 100 psia, a terminal impedance, 

Z t = 3.91 X 10 5 lb-sec/ft s . The results shown in Figures 78 and 79 reveal 
that for no external structural acceleration inputs, the feedline with the 
two 90° elbows behaves as a simple, straight feedline of the same total 
length. 

Line with Bends and Structural Excitation 

Figure 80 presents the feedline geometry and method of excitation 
used to experimentally determine the effects of structural-hydraulic 
coupling. The objective of this test was to' verify the analytical model 
described in Section V. 4 for a feedline with external structural accelera- 
tion applied at discrete locations. The mounting stiffness of the simu- 
lated support structure was quite high; i. e. , the mounting structure was 
almost perfectly rigid. The spring rate, k, of the mounting structure 
was determined from the elastic properties of the stainless steel shaft 
of the hydraulic actuator which was attached directly to the feedline 3 
inches from the second elbow, as shown in Figure 80. The structural 
damping of the mounting structure (shaft) was.zero for all practical pur- 
poses, and was thus programmed into the computer simulation. For the 
frequency range tested, the stiffness term, k, dominated such that the 
terms defined by Equation (114) became 

a/ => (Z J s) sinh yL 



As Z e sinh yL 
k 


a" => A/s (1 - cosh YL) 



A 3 s (1 - cosh YL) 
k ’ 


(155) 


and the relation (Equation 112) between the applied structural acceleration 
to the actual line acceleration, a^(s), simply became 
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Figure 79. Phase Angie Between the. Terminal Pressure and 
Pulsar Flow Perturbations for an Aluminum 
Line with Two 90° Bends 
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Figure 80. Feedline Configuration for Structu ral~Hycf raulic Coupling 
Tests 
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and the problem was essentially reduced to that of an accelerated line, 
•given- by Equation- (-1-08) * Although the- stiffness -term-dominated, -the ex- 
periment was modeled analytically using the more complex mounting 
stiffness model. Equation (113). Close examination of subroutine eight 
(given in Appendix D) reveals that the general computation method requires 
two horizontal line segments attached to the vertical line segment, one at 
each end. The problem under consideration, however, had only one hori- 
zontal segment. To properly program this problem, an additional hori- 
zontal line element was assumed to exist between the end of the vertical 
segment and the terminal resistance element. This extra horizontal seg- 
ment was programmed to have a length equal to zero. The addition of 
this imaginary element was necessary for subroutine eight to properly 
function. 

The initial data collected for this test configuration bore no re- 
semblance to the theoretically predicted response, and considerable time 
was taken to carefully re-examine the analytical model and the computer 
code. Re- consideration of the test procedures revealed that the actual 
input accelerations measured during the test were not the only source of 
structural-hydraulic excitation. Prior to the experiment it was assumed 
that the only excitation applied to the system would be that input created 
by the hydraulic actuator and would essentially invoke a vertical excita- 
tion of line segment 3-4 and bending of the horizontal line segment 2-3. 
Further tests revealed, however, that there was actually a dependent 
horizontal acceleration of considerable magnitude of the horizontal line 
segment; this acceleration was actually measured with an accelerometer 
mounted at the first elbow (point 2) and at certain frequencies was com- 
parable in amplitude to the measured vertical input excitation applied by 
the actuator. The experimental data were therefore the result of two 
simultaneous excitations of the feedline. To circumvent this problem, 
the feedline was attached rigidly to the support tower structure just up- 
stream of the first elbow. This bracing was very rigid, and prevented 
the secondary horizontal excitation when the second test was performed. 

The actual structural input was now very near that measured in the 
vertical plane produced by the hydraulic actuator. 

The terminal pressure/input acceleration is presented in Figure 81 
for a feedline with two bends and excited through structural-hydraulic 
coupling. The lengths of the line segments are given, in Figure 80; the 
internal fluid was water at a nominal line pressure of 60 psig at the line 
terminal. The terminal impedance was produced by the 7/l6-inch dia- 
meter orifice; the steady-state mean pressure drop was 60 psid producing 
a linearized terminal resistance, Zt = 2. 95 X 10 5 lb- sec /ft 5 . The excitation 
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Figure 81. Frequency Response for a Line with Bends and with Structural 
Excitation 
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amplitude was measured with the velocity transducer attached to the 
actuator shaft; the experimental velocity amplitudes were transformed 
to acceleration amplitudes by 

| ajg(tM) | = 0) | V £ (w)| (157) 

Figure 81 indicates that the analytical model does predict the 
amplitudes reasonably well past the first anti-resonance, as well as pre- 
dict the frequency at which this anti- resonance occurs. The scatter in 
amplitude measurements out to approximately 20 Hz is attributed to the 
low level pressure signals produced on the oscillograph trace, even 
though the system instrumentation was operated at maximum sensitivity. 
As yet we have been unable to explain the discrepancy between theory and 
experiment over the range from 40 Hz to 65 Hz, The analytical model 
did predict the resonant peak near 68 Hz; the experimental data show this 
peak to be at 70 Hz, The amplitude at resonance is, however, less than 
that predicted by theory. Between 70 and 82 Hz the system behaved very 
similar to the response predicted by the analytical model. The measured 
phase relations between the terminal pressure and acceleration inputs 
are presented in Figure 82, and the same general comments expressed 
above regarding the amplitude versus frequency response can be applied 
here also. 
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FREQUENCY, HERTZ 

Figure 82. Phase Angle vs Frequency for a Line with Bends and with 
Structural Excitation 
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VIII. CONCLUSIONS AND RECOMMENDATIONS 


VXII.1 -Conclusion s. 


The objective of this study has been to develop an analytical model 
and corresponding computer program to allow for the study of disturbances 
of liquid propellants in typical engine feedline systems. This model has 
been completed and includes (a) the effect of steady turbulent mean flow, 

(b) the influence of distributed compliances, such as dissolved ullage 
gases and flexible walls, (c) the effects of local compliances, such as 
cavitation regions and complex side branches, and (d) various factors 
causing structural-hydraulic coupling, such as bends, mounting stiffness, 
forces changes in line length and bellows or PVC joints. The computer 
program has been set up such that the amplitude and phase of the terminal 
pres sure /input excitation is calculated over any desired frequency range 
for an arbitrary assembly of various feedline components. A user's 
manual has been prepared and is attached as Appendix D of this Final 
Report. 


In addition to the development of the generalized computer code for 
the analysis of liquid rocket propellant feedline dynamics, an experimental 
test program has been conducted to verify many of the assumptions and 
models used in the computer simulation. Specifically, tests have been 
conducted to determine the validity of the analytical techniques used to 
model the effects of 

(a) steady, turbulent mean flow 

(b) linear and nonlinear resistance elements 

(c) blocked lines 

(d) elastic walls 

(e) dissolved gases 

(f) complex side branches 

(g) local compliances 

(h) bends-, and 

(i) structural-hydraulic coupling with mounting stiffness. 

In summary, the following conclusions can be drawn from the experience 
gained from both the analytical and experimental phases of this program: 
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(1) Investigation has shown that the predominant effect of turbu- 
lence is to increase the spatial attenuation at low frequencies; at high 
frequencies the laminar and turbulent attenuations coincide. The effect 
of turbulence also has a tendency to reduce the phase velocity at low 
frequencies; however, this effect can be neglected for all feedlines where 
the parameter Wr|/v > 10 4 , which includes virtually all cases of interest. 
An additional factor, r c > has been added to the laminar propagation 
operator to account for the turbulent attenuation contribution. 

(2) The terminal pressure amplitudes at resonance are affected 
more by the value of the line terminal impedance (resistance) than any 
other single factor. 

(3) The use of a linearized model for nonlinear type resistance 
elements provides satisfactory results for cases where the mean flow 
rate, Q, is considerably greater than the dynamic flow perturbation 
amplitude, Qd . 

(4) It has been concluded that the effect of the wall compliance 
(or elasticity) can be correctly modeled by the classic Kortew.eg correc- 
tion to the phase velocity; the attenuation for the elastic wall is virtually 
indistinguishable from the attenuation of a rigid wall. We have also con- 
cluded that for typical feedline problems, the axial wall stiffness effect 
will not permit real wave propagation of the higher order modes since 
the amount of energy which is fed into the higher order modes at the line 
termination will not permit sustained coupling of the fluid and the wall 
for these modes. 

(5) For the homogeneous distribution of very small gas bubbles 
entrained in the liquid propellant, the net effect over the frequency range 
of interest is to lower the phase velocity of the mixture. Experimental 
data indicate that there is also greater damping with increasing mass 
ratio than predicted by theory. 

(6) Small diameter gas bubbles entrained in the liquid reduce the 
phase velocity in the mixture by the amount predicted by an isothermal 
model for the gas speed of sound rather than the adiabatic model. 

(7) The model for a complex side branch such as an accumulator, 
predicts response amplitudes somewhat greater than those obtained ex- 
perimentally. The experimental resonant frequencies were shifted in 
accordance with theory, but agreement with theory becomes poorer as 
the terminal gas volume is increased. 
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(8) The analytical model used for a local compliance (large local 
gas volume) agrees very well with the data from the experimental' tests 
conducted. The experimental data also indicate that for cases where the 
gas bubble resonant frequency is far above that encountered through 
excitation, a more simple model for the compliance would be satisfactory; 
that is, a model which neglects the bubble dynamics. 

(9) A line with bends or elbows behaves exactly like a straight 
line of equal total length when the method of excitation does not involve 
structural acceleration of the line. 

(10) The model for the structural-hydraulic coupling of a feedline 
with bends (which provide local impedances) adequately predicts the line 
resonant frequencies, although the amplitudes measured were considerably 
different from the theoretical values. However, the specific example 
tested produced poor agreement with theory over a frequency range of 
40 to 65 Hz, and no satisfactory explanation for this behavior has been 
found. 


(11) The initial tests conducted involving structural- hydraulic 
coupling demonstrated the need for a thorough knowledge of the exact 
modes of excitation applied to a given feedline geometry. It has been 
found that certain excitations can produce dependent or secondary modes 
of excitation which must be properly treated before a correct analytical 
simulation is possible. 

VIII.2 Recommendations 


Based on the experience gained from this study of liquid rocket 
propellant feedline dynamics, recommendations for future work include: 

(1) The development of an analytical model and computerized 
subroutine to compute the speed of sound or phase velocity in two -phase, 
single- component fluids. The equilibrium model described in Section 
HI. 1 of this report would offer an excellent starting point, but would in- 
volve programming a vast amount of thermodynamic data for each liquid 
propellant under consideration. 

(2) Improvement of the model used to compute the bubble dynamics 
of large cavitation bubbles, where the bubble stiffness (or compliance) is 
also dependent upon the vaporization and condensation rates as the bubble 
undergoes successive pressure perturbations. 

(3) An analytical model should be developed to predict the effects 
of nonlinear terminal resistances for configurations or operational condi- 
tions where Q » Q d . The nonlinear distortion created at these conditions 
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could prove significant; this analytical model should also be examined 
experimentally. 

(4) Further experimental tests should be conducted to determine 
the validity of the model used to simulate structural-hydraulic coupling 
and the effects of mounting stiffness. 

(5) An analytical and experimental study should be conducted to 
study the apparent added damping due to the presence of entrained gas 
bubbles. The data obtained in this study indicate that the increased 
damping is significant and should not be neglected. 

(6) General modification of the feedline computer code described 
herein to allow more versatility in the use of this program. For example, 
the present program can currently handle only one mode of excitation at a 
time, whereas, in reality, there can be numerous sources of simultaneous 
feedline excitation. 
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APPENDIX A 

RATIONAL APPROXIMATE MODEL FOR 
DISTRIBUTED PARAMETER SYSTEMS 


A- 1 



The four -terminal representation of a transmission line that 
has been utilized throughout this report is an exact model for the 
zeroeth mode transfer function equations. The value of this repre- 
sentation for frequency response analyses has been proven. In 
many problems, -however, -the time domain-solution is required; - and" 
the Laplace transform representation constitutes a convenient tech- 
nique for obtaining the space- dependent solution independently of the 
time -dependency. In the final analysis, the Laplace inversion inte- 
gral must be evaluated. When hyperbolic functions are involved in 
the inversion, as they are in the distributed parameter line model, 
an infinite series of time -dependent terms is obtained due to the un- 
bounded number of zeros of the hyperbolic functions. In many cases, 
closure of this series is extremely difficult, and the desired accur- 
acy does not warrant laborious mathematical operations. The ob- 
vious conclusion is that there is a definite need for a simple, accu- 
rate, approximate technique for obtaining the time domain solution 
from the frequency domain representation. 

The rudiments of such a technique were first proposed by 
Oldenberger and Goods on^, and were later refined and put in a 
practical form by Gerlach^' This technique is based upon ex- 
panding the hyperbolic functions into an infinite product of second- 
order polynomials, i. e. , 

courts) = n n Q { 1+ ^ + ^} (A.i) 

00 2 

sinhr(s) =r(s) n {i+^j^ + .-S-} (a.2) 

The coefficients, £ cn , U) 0n , £ sn , w sn are evaluated by solving for the 
value of Laplace operator, s„ at the zeroes of the hyperbolic func- 
tions. These quantities are conveniently catalogued in Figures A-l 
to A-4. To evaluate any one of these coefficients, it is sufficient to 
calculate the damping number of the line, D n . The damping number 
and the desired term number in the expansion completely define the 
four coefficients. It may be shown that for approximation purposes 
the term F(s), in the expansion for the hyperbolic sine, can be 
represented adequately by sL/c 0 . Characteristic impedance can 
be approximated by p 0 c 0 . 

With these expansions, it is now possible to approximate the 
hyperbolic functions to any desired degree of accuracy with expres- 
sions that can easily be (1) inverted into a time solution using 
standard transform pairs or (2) integrated in the time domain to 
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Figure A-l. Variation of the Approximate Model Parameter 
F cn With Axial Damping Number 
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Figure A- 2, Variation of the Approximate Model Parameter 
C C n With Axial Damping Number 
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Figure A- 3. Variation of the Approximate Model Parameter 
F sn With Axial Damping Number 






obtain the real time solution. For example, consider a transfer 
function in which a one-term approximation for the hyperbolic 
functions produces 


P(s) 

P 0 C 0 V 0 



l+‘ 


176 

66 


s + 



Standard inversion pair yields 


P(t) 

Po C o Vo 


1. 5e 


-5. 8t 


sin (66t) 


{A. 3) 


(A. 4) 


As a second example, consider a transfer function of the type 


P(s) = -Z c (s) tanh T{s) V(s) (A. 5) 

A one-term approximation in the model yields 



interpreting the Laplace operator as a time derivative, the follow- 
ing equation is obtained: 



CQ 


Wo 


jd_ JL cf_\ 

dt + W co 3 dt? J 


p(t) 


-P 0 L 


dv(t) 

dt 


(A. 7) 


With a second equation relating pressure and flow velocity, a stan- 
dard numerical integration easily can be obtained on the computer. 


The accuracy of this model is summarized below: 

(1) One term of the model well approximates the hyper- 
bolic operators up to the first critical frequency. 

(2) Two terms improve the approximation up to the first 
critical point and roughly (not well) approximate the 
hyperbolic operators up beyond the second critical 
frequency. The use of more terms would improve the 
results near the second critical frequency. 

(3) A one-term approximation for the hyperbolic function 
is generally adequate for most engineering problems. 
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APPENDIX B 

SPEED OF SOUND IN A LIQUID CONTAINING 
A HOMOGENEOUS DISTRIBUTION OF BUBBLES 


B- 1 



Single- Component, Two-Phase Mixture in Equilibrium 


The speed of sound in a pure, two-phase substance can be 
expressed by the standard form 



(B.l) 


or, in terms of specific volume rather than density. 


a 3 = 



(B.2) 


Following the approach of Gouse and Brown®, it can be shown that 
the partial derivative in Equation (B.2) can be expanded such that 


a 2 = 



(B.3) 


Since constant pressure and temperature lines are coincident in the 
two-phase region beneath the vapor dome. 



or in finite difference form 

/ 5 s\ _ s e - s t 

\ a v/ p Vg-v f 


(B.4) 


(B.5) • 


The subscripts, f and g, refer to the saturated liquid and vapor 
state, respectively. Substituting Equations (B.4) and (B.5) into 
Equation (B.6) 



Equation (B.6) can be developed into a form that is suitable 
for use with a standard temperature- entropy chart. The quality of 
a two-phase state, X, is defined as 

X - M. V /(M X+ M V ) (B.7) 

and the specific volume for such a state can be written in terms of 
quality as 
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V = Vf + x(v g - v f ) 


(B.8) 


Substituting this expression into Equation (B.6), the speed of sound 
can be written as a function of quality 





stf 



(B.9) 


In finite difference form, the speed of sound in a single- component, 
two-phase mixture in thermodynamic equilibrium becomes 


a 




(B.10) 


Single- Component or Two- Component, Two-Phase Mixture With 
Constant Quality 

The above section described the acoustic velocity in an 
equilibrium mixture of a liquid and its vapor phase. In this section, 
the speed of sound is analyzed for a constant quality mixture of a 
liquid and a distinctly different gas, but the results also apply to 
single -component mixtures as long as the vaporization and conden- 
sation rates are low enough that no phase change occurs as a result 
of pressure disturbances {constant quality is assumed at all times). 
In the development of a constant quality model for the acoustic velo- 
city in a liquid-gas mixture, the following assumptions are made: 

(1) The mixture is homogeneous in phase composition. 

(2) The mass of each phase remains constant. 

(3) The gas behaves as a perfect gas with the appropriate 
equation of state and constant specific heats. 

(4) The liquid compressibility is ignored. 

(5) The gas and liquid phases are always at the same 
temperature. 


The subsequent derivation is based on total volume as opposed to the 
specific volume approach that was utilized in the previous section. 
Assuming an adiabatic process. 


c 


2 



or 



(B.ll) 
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For a two- component mixture. 


P = 


M g 

v J e+v g 


_ _a 

c 3 " dp \ V£+ V g ) 


^ = [(Va+V 6 ) ^ (M^+Mg) - 

J_ (M a± M A 1 .. . , _1 

c 2 " ' Wi+v E ) Yjg+y, dp 


(Mi+M g )^ (V^+ V g )]/ (Yjj+V g ) 2 
(V^+V g ) 



(B.12) 


(B.13) 


From the definition for density, 

V = M/p 


_ Ml , Ml _ Pe + M g PA 
PA P 6 PAP e 


(B. 14) 
(B. 15) 


The mass ratio, $, is now defined as the mass of vapor/mass of 
liquid. 


$ = M s /Mi 

Substituting the mass ratio into Equation (B.15), 

v ,+ v = .¥ 4 .(£i±iM 

1 6 P 6 Pi 


(B.16) 


(B.17) 
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Thus, the speed of sound in a constant quality mixture becomes- 




— P B 
P g +4Pi 


_ 1 _ _ 1 _ 

cf 


. _i£i_ 

px+§p g 


J. 

Pg 



(B. 18) 


As mentioned in Section III, the isothermal velocity of sound 
might be more appropriate for two - component, ' two- phase mixtures 
rather than the adiabatic (isentropic) velocity of sound* 

Pies set and Hsieh^ have shown that bubble dynamics are 
governed by an isothermal process at low excitation frequencies and 
an adiabatic process at higher frequencies, providing there is a uni- 
form temperature distribution within the gas bubbles. However, 
this sharp division of thermodynamic behavior is not necessarily 
the case because the heat conduction rate of the gas is considerably 
smaller than that in the liquid. The liquid has a large specific heat 
and thermal conductivity, and behaves as a heat reservoir. Conse- 
quently, in the liquid adjacent to the gas -liquid interface, there are 
no changes in the gas temperature during compression. In the 
center of the bubble away from a substance having a high specific 
heat, the gas follows the adiabatic equation of state. Therefore, the 
overall bubble follows a polytropic process. However, for very 
small bubbles, the heat diffusion into the liquid is so rapid that 
temperature changes in the bubble can not take place, and the pro- 
cess is essentially isothermal. 


For this reason, the propagation of sound in a liquid with a 
homogeneous distribution of bubble s _ {b~a sic assumption No. 1) indi- 
cates a speed in agreement with isothermal speed of sound. To 
compute the isothermal speed of sound in the mixture, the acoustic 
velocity of the gas in Equation (B.13) should be changed to the iso- 
thermal value, or 


c g =Vgo YRT 1 

where Y = 1. 0 rather than 1.4. 


(B.19) 
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APPENDIX C 
PARALLEL LINES 


C-l 



Consider the assemblage of parallel lines shown in Figure 
C-l, Although this configuration is not used frequently in pro- 
pellant feed systems, the pressure-flow relationships for this com- 
ponent have been developed and included in the computer code for 
application in special situations. Assuming that flow losses due to 
elbow- resistance- a-r-e negligible^ a- -standard- -four -terminal- pressure- 
flow equation can be written for the i-th line. For this particular 
problem, however, it is more convenient to express that pressure- 
flow equation as 


Qi t = Zcj 1 At (Pi coth I* - PsCschTi) 

Q 3l = Zcj 1 Ai (P x csch Ti - P 2 cothri) 


(C.l) 


The terminal pressures, and P 2 , are common to all line seg- 

ments* The continuity equation requires that 


n n 

E Qii = Qi and E Qs i = Qs 

i=l i=l 


(C.2) 


Summing Equations (C.l) and utilizing Equations (C.2) yields 


n n 

Ql = Px Y Z Cl X Ai coth r t - P 2 £ Zoj 1 Ai csch r t 
i=l i=l 


-l 


n n 

Q 2 = Pj. E Z Ci X Ai csch Ti - P 2 ^ Z Ci ‘Ai coth II 
1=1 i=l 


(C.3) 


Equations (C.3) can be easily rearranged into the standard form 
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where 
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Figure C-l. Parallel Line Model 
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The line length that appears in the expression for the propagation 
operator is taken to be the developed length measured along the 
pipe centerline. 
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APPENDIX D 

USER'S MANUAL FOR GENERALIZED 
FEEDLINE PROGRAM 
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D. 1 Program Organization 


This manual describes the mechanics of the Generalized 
Feedline Computer Code. It is assumed that the user has had pre- 
vious experience in FORTRAN IV coding. 

The FORTRAN IV source deck is composed of a main or 
controller program, individual subroutines for each line component 
and several special purpose subroutines which either perform 
specific mathematical operations or. serve as function evaluators. 
The subroutine approach was chosen because of its distinct advan- 
tages in the areas of program debugging, modification and expan- 
sion. Using the Control Data Corporation 6400 system, the source 
deck is compiled each time the program is submitted to the terminal 
for execution, thus eliminating tape handling. The function of each 
routine is presented below in the order in which it will appear in the 
FORTRAN listing which is included in the next section. 

Controller: 


The controller program performs the following functions: 


(1) Read-in of all pertinent data for a particular line 
configuration; 

(2) Execute the calling of subroutines in the proper sequence; 

(3) Evaluate the magnitude of the transfer function at each 
frequency; 

(4) Increment frequency, the independent variable, within 
a prescribed range; 

(5) Perform all output operations. 


Subroutine ONE : Calculates matrix elements for a straight duct 

of finite length and constant cross-sectional 
area having no external excitation. 


Subroutine TWO: Calculates matrix elements for a single cavi- ' 

tation bubble including the effects of radiation, 
thermal and viscous damping. 

Subroutine THREE: Calculates the matrix elements of a pressure- 

volume compensator (PVC) joint including the 
effects of internal friction. 
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Subroutine FOUR; 
Subroutine FIVE: 

Subroutine SIX : 

Subroutine SEVEN: 

Subroutine EIGHT : 

Subroutine NINE: 

Subroutine TEN : 

Subroutine ELEVEN : 

Subroutine TWELVE 
to FIFTEEN : 

Subroutine SPEED: 


Subroutine DMAT: 
Sub r outine BMAT: 


•Establishes the continuity requirements for 
a side branch pulser. 

Calculates the matrix elements for a straight 
line which undergoes rigid body motion as a 
result of an external velocity excitation. 

Calculates the matrix elements for several 
lines in parallel all of which have a common 
entrance -exit point. 

Calculates the matrix elements for a bellows 
with gas trap liner. Relative motion between 
the ends of the bellows is permitted. 

Calculates the matrix elements for a line with 
mounting stiffness having a structural accelera- 
tion applied to the viscoelastic support. 

Calculates the matrix elements for a line 
undergoing forced changes in length. 

Calculates matrix elements for a line with 
mounting stiffness using the impedance method. 

Calculates the matrix elements for a com- 
plex side branch with laminar flow. 


Available for future expansion. 

Calculates the speed of sound in a pure liquid 
or a liquid with a homogeneous distribution of 
an ullage gas. A mixture of a liquid and its 
own vapor may also be considered. However, 
in the latter case, the results are not precise 
because the theoretical model does not in- 
clude the thermodynamic effect of phase change 
at the bubble surface. Finally, the speed of 
sound may be corrected for wall elasticity. 

Calculates the resultant matrix D in 
Equation (122). 

Calculates the resultant matrix B in 
Equation (123). 
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Subroutine J1ORJ10 : - Evaluates Bessel functions Jj.(Z) and J 0 (Z) 

for complex Z. Series expansion is used 
for moderate values of the argument; an 
asymptotic approximation is used for large 
arguments. 

D. 2 Program Listing 

This section contains a listing of the program source deck. 
Program control cards, which are described in Section D. 4, are 
applicable to CDC 6000 series computer systems. 

D. 3 Program Input Package 

Due to the generalized nature of this computer program, the 
length of the input package is variable and is directly related to the 
number of components in the feedline model. The first card in the 
package is always an alphanumeric header card on which the user 
can punch, pertinent information such as run number, configuration 
number, etc. , to be printed as the first line of program output. 
Fifty-five (55) Holerith fields have been allocated for the informa- 
tion on the header card. The actual data input is a combination of 
integer or floating point information. Integer data conforms to a 
(2413) format, while the floating point data is read in with a (6E12.6) 
format. The following data must be supplied in the order and with 
the units shown: 


CARD 1 Header card 


CARD 2 


NELM, dimensionless Designates the number of 

components in the line. 

JBNUM, dimensionless Number of matrices that 

must be summed to arrive 
at matrix B, Equation (123). 


NGAS, dimensionless Indicates presence (NGAS = 1) 

or absence (NGAS = 0) of dis- 
solved gases in propellant. 


NELAST, dimensionless Includes (NELAST - 1) or 

excludes (NELAST = 0) line 
elasticity effects in speed 
of sound calculation. 

NEXCITE, dimensionless Element type producing the 

primary line excitation. 
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PROGRAM CONTROL (INPUT , OUTPUT, TAPE bO=INPUT) 

NELMsNUMBER OF DISCRETE ELEMENTS IN THE FEEDLINE, I. E. , LINES, 
BELLOWS, PVC JOINTS, ETC. 

NELM INTEGER CONSTANTS ARE READ IN AND DESCRIBE THE TYPE OF 
ELEMENTS IN THE ORDER THEY APPEAR IN THE LINE BEGINNING AT THE 
OUTPUT OF THE FUEL TANK AND TERMINATING WITH THE COMBINED INPUT 
IMPEDANCE TO THE TURBOPUMP, IN JECTOR AND ENGINE. FOR EXAMPLE, 
ITYPE( J)=M, JsELEMENT NUMBER, M=ELEMENT TYPE 
M=1 SIMPLE LINE WITH NO EXTERNAL EXCITATION 
M=2 CAVITATION BUBBLE 
M=3 PRESSURE VOLUME COMPENSATOR 
Ms* SIDE BRANCH PULSER 

M=S RIGID BODY MOTION OF A SIMPLE LINE, VELOCITY EXCITATION 
M=b NPAR LINES IN PARALLEL WITH A COMMON INPUT-OUTPUT POINT 
Ms? BELLOWS WITH RELATIVE MOTION BETWEEN THE ENDS 
M=S LINE WITH MOUNTING STIFFNESSCOPTION 1) 

M=q FORCED CHANGES IN LINE LENGTH 

M=10 LINE WITH MOUTING STIFFNESSCOPTION 2) 

M=ll COMPLEX SIDE BRANCH 

M=12 

M=13 

M=If 

M-15 

JBNUM=NUMBER OF MATRICES THAT MUST BE SUMMED TO CONSTRUCT MATRIX B 

B=B(1)+B(2)+ +BCJBNUM) 

JTERM(J) IS THE NUMBER OF SUBMATRICES IN B(J) 

K( J,M) IS AN ARRAY THAT DESCRIBES THE TYPE OF MATRICES THAT ARE TO 
BE MULTIPLIED TOGETHER TO CONSTRUCT THE COMPONENTS OF MATRIX B. 

M TAKES ON VALUES FROM 1 THROUGH JTERM(J) 

FOR EXAMPLE, IF B(2) =DS*D<**C3 THEN K(2,l}s5,K(2,2)=<t,K(2,3)=3 
NGAS=Q, NO DISSOLVED GASES IN THE PROPELLANT 

NGASsi, DISSOLVED GASES IN PROPELLANT WITH GAS-TO-LIQUID MASS 
FRACTION, PHI 

NELASTsO, SPEED OF SOUND IN PROPELLANT IS BASED ON AN INFINITELY 
RIGID TRANSMISSION LINE 

NELAST=1, SPEED OF SOUND IN PROPELLANT REFLECTS LINE ELASTICITY 
(WALL MODULUS EWALL AND WALL THICKNESS HWALL) 

NEXCITE=NELM PRODUCING THE PRIMARY LINE EXCITATION 
BRCOMP=COMPLIANCE OF SIDE BRANCH ELEMENT 
8RL=LENGTH OF SIDE BRANCH ELEMENT IN FEET 
BDIAM=DIAMETER OF SIDE BRANCH ELEMENT IN INCHES 
BRRsRESIST ANCE OF SIDE BRANCH ELEMENT 
********************************* 
COMPLEX ET A, ENUM, DENOM, TRANS, B,D, XI, ARG,RJ, GAMMA, ZC,COSHG,SINHG,DD 
1, CC, TCOTH, TCSCH, BB, Z1,Z2, XI* X2,X3,Xf , QUADRAT, ALPHAP, BETAP, ALPHADP, 
2BETADP,GOFS,SLCW,COSHCW,SINHCW,COEFF,ALPHA1, ALPHAS 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION I TYPE (15), JTERMC15) ,K(1S,1S) ,EL(1S) ,RADIUS(1S) ,RBUB(1S) , 
1FPVC(15),AKBPVC(15),AKCPVC(15),PARLEN(1S,5),PARRAD(15,5),FREQ(200) 
DIMENSION B(2,1),DC2,2),SIZE(200),SIZEDB(200),ANGLE(200),FBEL(15), 
1C0MPLY(15),AKBEL(15) ,BSIGN( 10), DAMPER CIS), SPRINGK (15) 

DIMENSION DD(2,2,15),CCC2,1,15),BB(2,1,15),RADSEC(200),AREA(15) 
COMMON PI,RADIUS,OMEGA,ENU,ETA,EL,CO,RHOO,THERMK,RBUB,'CPCV,PO,VISC 
1,FPVC, AKBPVC, AKCPVC, NPAR,PARRAD,PARLEN,BULKMOD,RHOLIQ,NGAS,GASMW 
COMMON GAMG AS, G, TO, PHI, COMPLY, NELAST, EWALL, HWALL, NELM, JBNUM, JTERM, 
1K,DD,CC,D,B,BB,CP, J, AREA, VMEAN,FBEL, AKBEL, BSIGN, DAMPER, SPRINGK 
COMMON G1,G2,RH0WALL,ZX,ZY 
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COMMON BRCOMP,BRL,BRDIAM 

C ********************************* 

1 READ 1000 

IFCEOF, b0)5, 10 
5 STOP 
10 CONTINUE 

READ 1005,NELM, JBNUM/NGAS, NELASTr NEXCITE 
READ 100S,CITYPE(J)f JslrNELM) 

READ 1 OJD 5 * C.JT.E R M.( JXr- J.= L,- J B N U M-) 

DO 15 J=1,JBNUM 
KOUNTsJTERM(J) 

15 READ 1005, (K C J , M) , M=1 , KOUNT) 

READ 1020,(BSIGN(J), J=1,JBNUM) 

READ 10BO,HERTZI,DELHZ1,F11IM,DELHZ2,HERTZF 
READ 1020, SIGN, TERMZ , RHOLI Q , THERMK , PO , TO 
READ 1020,EWALL,HWALL,ZINPUT,CPCV,BULKM0D,PHI 
READ 1020, VISC, GASMW,GAMGAS, CP, VMEAN 
C LOOP FOR READ-IN OF LINE ELEMENT DATA 

PI=3.1tl5S2? 

DO 55 J=1,NELM- 
INDEX=ITYPE(J) 

GO TO (20,2S,30,3S,t0,t5,50,55,b0,b5,70,75,80,85,S0)»INDEX 
20 READ 1020,EL(J),RADIUS(J) 

AREA(J)=PI*RADIUS(J)**2/ltt. 

AREA1=AREA(J) 

RAD1=RADIUS(J)/12. 

GO TO 55 

25 READ 1020,RBUB(J) ,RADIUS(J) 

GO TO 55 

30 READ 1020, FPVC ( J) , AKBPVC ( J) , AKCPVC ( J) 

GO TO 55 
35 GOTO 55 

tO READ 1020,EL(J),RADIUSCJ) 

AREA(J)=PI*RADIUSCJ)**2/ltt. 

GO TO 55 

t 5 READ 1005, NPAR 

READ 1020 ,(PARLEN(J,n,PARRADCJ,I),I=l,NPAR) 

GO TO 55 

50 READ 1020,FBELCJ),C0MPLYCJ),AKBEL(J) 

GO TO 55 

55 READ 1020,£L(J),RADIUS(J),DAMPER(J),SPRINGK(J) 
AREA(J)=PI*RADIU$(J)**2/ltt. 

GO TO 55 

bO READ 1020,EL(J),RADIUS(J),G1,G2,RH0WALL 
AREA(J)sPI*RADIUS(J)**2/ltt. 

GO TO 55 

b5 READ 1020,EL(J),RADIUS(J),ZX,ZY 
AREA(J)=PI*RADIUS(J)**2/l*t. 

GO TO 55 

70 READ 1020,BRL,BRDIAM,BRC0MP 
GO TO 55 
75 READ 1000 
GO TO 55 
80 READ 1000 
GO TO 55 
85 READ 1000 
GO TO 55 
50 READ 1000 
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«15 CONTINUE 

C ********************************* 

G=32.17f0<*q 
ETA=(0.,1.) 

IFREQsl 

FREG(IFREG)=HERTZI 

C CALCULATION OF MATRIX ELEMENTS 

55 OMEGA=2.*PI*FREQ(IFREG) 

DO 175 J=1 » NELM 
INDEX=ITYPE(J) 

GO TO (100,105,110,115,180,125,130,135,1*0,1*5,150,155,150, 155,170 
11, INDEX 
100 CALL ONE 
GO TO 175 
105 CALL TWO 
GO TO 175 
110 CALL THREE 
GO TO 175 
115 CALL FOUR 
GO TO 175 
120 CALL FIVE 
GO TO 175 
12S CALL SIX 
GO TO 175 
130 CALL SEVEN 
GO TO 175 
135 CALL EIGHT 
GO TO 175 
1*0 CALL NINE 
GO TO 175 
1*5 CALL TEN 
GO TO 175 
150 CALL ELEVEN 
GO TO 175 
155 CALL TWELVE 
GO TO 175 
IbO CALL THIRTEN 
GO TO 175 
155 CALL FOURTEN 
GO TO 175 
170 CALL FIVETEN 
175 CONTINUE 

C ********************************* 

CALL DMAT 
CALL BMAT 

ENUM=Ba,l)*CDC2,l)*ZlNPUT+0C2,2))-B(2,l)*(0(l,l)*ZINPUT+DCl,2)) 
DEN0M=DC2,i)*ZINPUT+D(2,2)-CD(l,l)*ZINPUT+D(l,2))/TERMZ 
TRANS=SIGN*ENUM/DENOM 
SIZE(IFR£Q)=CABS (TRANS) 

IFCNEXCITE.EQ.f )G0 TO l?b 
IF(NEXCITE.EQ.5)G0 TO 177 
IF(NEXCITE.EQ.7)G0 TO 17? 

IF(NEXCITE.EG.q)GO TO 177 
IF(NEXCITE.EQ.B)GO TO 178 

l?b SIZE0B(IFRE0)=20.*AL0G10CSIZECIFREQ)*AREA1/SQRT(RH0LIQ*BULKM0D)> 

GO TO 175 

17? SIZEDB(IFREQ)=20.*ALOG10(SIZE(IFREG)/SGRT(RH0LIQ*BULKM0D)) 

GO TO 175 
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178 SIZEDB(IFREQ)=B0.*AL0G10(SIZE(IFREQ)/(RH0LIQ*RAD1)) 

179 CONTINUE 

ANGLE(IFREQ)=ATAN2(A1MAG( TRANS) , REAL (TRANS) )*57. 29578 

180 ANGLE (IFREQ)=ANGLE(IFREQ) -180. 

185 RADSEC (,IFREQ)=OMEGA 

IF(FREQUFREG).LT.HERTZF) GO TO 190 
GO TO 80S 
190 TFREQ=IFREQ+1 

IF(FREG(IF.REQ-1)-F1LIM)19S,19S,200 
1 9'S' F 'R'E'Q ( I F R E Q ) s F R E Q ( I F R E Q - i ) + D E L H Z 1 
IFINAL=IFREQ 
GO TO 99 

200 FREQ (IFREQ)=FREQ( IFREQ-D + DELHZ2 
IFINAL=IFREQ 
GO TO 99 
205 PRINT 1000 
PRINT 1025 
PRINT 1030 
00 250 1 = 1 , IF INAL 

250 PRINT 1090,RADSEC(I),FREQ(I),SIZE(I),SIZEDB(I),ANGLECn 
GO TO 1 

c ********************************* 

1000 F0RMATCSSH1 ) 

1005 FORMAT (2913) 

1020 FORMAT (bEl2.b) 

1025 FORMAT (1H0,39X,15HN0N-DI MENS I0NAL,5X,5HPHASE) 

1030 FORMAT (19H OMEG A-R AO/SEC , 4X , 7HFREQ-HZ, fX , 9HTRANS-P/K , feX, 8HTRANS-DB 
1, bX, 9 H ANGLE- DEG) 

1090 F0RMAT(1X,F9.1,F15.1,F13.3,F19.2,F15.1) 

1050 F0RMATClHl,Elf.b,bElS.b) 

END 
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SUBROUTINE ONE 

THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE TRANSFER MATRIX FOR 
A STRAIGHT LINE OF LENGTH L WITHOUT EXTERNAL EXCITATION 
COMPLEX ETA f ENUM,DENOM,TRANSr6,D,XI, ARG,RJ f GAMMA, ZC, COSHG, SINHG, DO 
lr CC,TC0TH,TCSCH,B6,Z1,Z2, XI , X2,X3,X«f , QUADRA T, ALPHAP, BET AP, ALPHADP, 
2BETAOPf GOFSf SLCW,COSHCWf SINHCWf COEFFf ALPHA1, ALPHA2 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION ITYPEC15), JTERM(iB), K( 15,15 ),EL(15), RADIUS (IS), RBUB (15), 
1 FPVCC 1 S), AKBPVCC15) , AKCPVCC15), PARLENClS, 5),PARRADas,5)f FREQ(eoO) 
DIMENSION B(2,1),D(2,2),SIZE(200),SIZEDB(200),ANGLE(200),FBEL(15), 
.1 COMPLY CIS ) , AKBELC1S) , BSlGN(lQ) , DAMPER(IS) , SPRINGK (15) 

DIMENSION DD(2,2,15),CC(2,l,15),BB(2,l,i5),RADSEC(2(lG),AREA(15) 
COMMON PI, RADIUS, OMEGA, ENUf ETA, EL, CO fRHOOfTHERMK, RBUBf CPC VfPQf VISC 
i,FPVC, AKBPVC, AKCPVC,NPAR,PARRAD,PARLEN,BULKMOD,RHOLIQ,NGAS,GASMW 
COMMON GAMGAS, G, TO, PHI, COMPLY, NELAST,EWALL,HWALL,NELM,JBNUM,JTERM, 
1K,DD,CC,D,B,BB,'CP, J, AREA, VME AN, FBEL, AKBEL , 8 S I GN , DAMPER, SPRINGK 
COMMON G1,G2,RH0WALL,ZX,ZY 
CALL SPEED 

REYN0LD=VMEAN*2.*RADIUS(J)/(ENU*12.) 

RT=lH‘t.*(2.*ENU*a.005 5*REYN0LD**n.85)/(RADIUS(J)**2) 

Zl=(OMEGA*EL( J)*ETA)/CO 
Z2=CSQRT((l.,0.)+R1/(OMEGA*ETA)) 

XI = SQRTCUMEGA/ENIJ)*CSQRT(-ETA) 

ARG=XI*RADIUS(J)/12. 

CALL J10RJU(ARG,RJ) 

GAMMA=ETA*OMEGA*EL(J)/(CO*CSQRT((1.,G.)-2.*RJ/ARG))+REAL(Z1*Z2) 

ZC = RHOn*Cn/(CSQRT( (.1 ., 0 .)- 2 .*RJ/ARG)) 
C0SHG=(CEXP(GAMMA)+CEXP(-GAMMA))/2. 

SINHG=(CEXP(GAMMA)-CEXP(-GAMMA) )/?. 

DDU , 1, J)sCOSHG 

00(1,2, J)=-ZC*SINHG/ AREA (J) 

DO (2, I, J)=-AREACJ)*SINHG/ZC 

DD(2,2, J)=COSHG 

RETURN 

ENO 
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SUBROUTINE TWO 

THIS SUBROUTINE CALCULATES THE COMPLIANCE OF AND TRANSFER MATRIX 
ACROSS A SINGLE BUBBLE 

CPC VsRATIO OF SPECIFIC HEATS FOR THE GAS IN THE BUBBLE 
CP=SPECIFIC HEAT AT CONSTANT PRESSURE FOR THE GAS IN THE BUBBLE 
COMPLEX ETA,ENUM, DENOM, TRANS, B , D , XI , ARG, R J, GAMMA, ZC/COSHG, SI NHG,DD 
1, CC,TCOTH,TCSCH,BB,Zl,Za, XI f X?/ X3,X*, QUADRAT, ALPHAP, 8ETAP, ALPHADP, 
aBETADP,GOFS, SLCW,CQSHCH,SlNHCW,COEFF, ALPHA! , ALPHAS 
COMPLEX ZSTRUCT,ALP.HA.,..B£XA 

DIMENSION ITYPEC15), JTERMU5), K C 15 , 15 ) , EL C 15 ) , R ADIUS CIS ) r RBUBClS), 
iFPVC(J5)f AKBPVC(I5)r AKCPVC(XS)fPARLENCl5r5)rPARRADCl5rS)rFREQ(200) 
DIMENSION B(2,l),D(2,23,SIZE(2D0),SIZEDBCa00),ANGLEC2GQ),FBELClS), 
J COMPLY CIS ),AKBELC IS), RSIGN( 10), DAMPER ( 15 ),SPRINGK (15} 

DIMENSION DDC2,8, 15),CC(E,1 ,15),BB(2»1, 15 ) , R ADSEC (200 ) , AREA (15) 
COMMON PI,RADIUS,OMEGA,ENU,ETA,ELfCO,RHOO,THERMK,RRUB,CPCV,PO,VISC 
l/FPVCr AKBPVC/AKCPVCf NRAR,PARRAD,PARLEN,BULKMOD,RHOLIQf NGAS/GASMW 
COMMON G AMGAS,G, TO, PHI, COMPLY, NEL AST, E WALL, HN ALL, NELM,JBNUH,JTERM, 
1K,DD,CC,D,B,BB,CP, J, AREA, VMEAN, FEEL, AKREL, 8S IGN, DAMPER, SPR I NGK 
COMMON G1,G2,RH0WALL,ZX,ZY 
GA5C0N=?7S.*CP*Cl.»l./CPCV) 

GASDENsl t l t f.*PO/(GAvSCON*G*(TO+‘lbO.}) 

DTHERMs fHERMK/CGASDEN*CP*G) 

PRO=SQRT(OMEGA*3feOO./(a.*DTHERM))*RBUB(J)/ia. 

IF ( PRO .GT. 35.)ao,in 

10 Tl = (SINh(2.*pRO) + SIN(a.*PRO))/(COSH(a.*PR[l)-COS(a.*PRa)) 
Ta=CSINH(2„*PRO)-SINC3.*PRa))/CCOSHCa.*PR0)-COS(a.*PRO)) 
T3=a.*PR0/(3.*(CPCV-l.)) 

Tf=(Tl-*l./PR0)/(T3 + ia) 

POLYsCPCV/CCl. + Ta/Ta^Q. + TM-^a) ) 

GO TO 30 
ao poly-cpcv 

30 BU8V0L = ^.*PI*RBUB(J)**3/(3.*1728.3 
CALL SPEED 

SPR ING-POL Y*P0*1*IH . / BUB VOL 
AMASSa=RHOO/ (f .*PI*RBUB( 

IF(PRQ»GT.35,}50, t H] 
fO BTHERM = T'f*SP RING/OMEGA 
GO TO 50 

50 8THERM=3.*(CPCV-l.)*SPRING/(a.*PRO*0MEGA) 

50 BRAD= AM ASSa*R8UB(J)* COMEG A**e)/(ie.*CG) 
8VIS=17?8„*VISC/(PI*R8UB(J)**3) 

pdamp=btherm+brad+bvis 

D0Cl,l,J}=Cl./0,3 

DDCl,2,J)s(o.,0.) 

00(2,1 , J)s-ETA*0MEGA/c-AMASS 2*0MEGA**2+ETA*BDAMP*0MEGA + SPRING) 
DDca,a, J)=(i.,n.) 

RETURN 

END 
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SUBROUTINE THREE 

THIS SUBROUTINE CALCULATES THE TRANSFER MATRIX FOR A PRESSURE- 
VOLUME COMPENSATORCPVC JOINT) 

COMPLEX ETA,ENUM,DENOM,TRANS,B,D,XI,ARG,RJf GAMMA, ZCrCOSHG, SI NHG,DD 
1,CC, TCOTH,TCSCH, BB,Z1, Z2, Xl ,‘X2,X3,Xt, QUADRAT , ALPHAP , BETAP, ALPHADP , 
2BETADP,GOFS,SLCW,COSHCW,SINHCW,COEFF, ALPHAl, ALPHAS 
COMPLEX ZSTRUCTr ALPHA, BETA 

DIMENSION ITYPE(15), JTERM (15 ) , K Cl 5 , 15 ) , EL( 15 ) , R AD IUS C 15) , RBUB ( 15 ) , 
1FPVC(15) ,AKBPVC(15) ,AKCPVC(15) , P ARLEN ( 15 , 5 ) , P ARR AD ( 1 5 , 5 ) , FREQ ( 200 ) 
DIMENSION B(2,1),D{2,2) ,SIZEC2G0) ,SIZEDB(200), ANGLEC200), FBELC15), 
1C0MPLYQ5) /AK6ELC15) ,RSIGN Qg) , DAMPER f 15) ,SPRINGKCl'5) 

DIMFNSION DD(2,2,1S),CC(2,1,15},BBC2, 1,15),PADSEC(200),AREA(15) 
COMMON PI , R ADIUS, OMEGA , EMU, FT A, EL, CO , RHOQ , THERMK , RBUB, CPCV , PO , V ISC 
1, FPVC, AKBPVC, AKCPVC, NPAR, PARR'AD, PARLEN, BULKMOD, RHOLIQ, NGAS, GASMW 
COMMON GAMG ASrG,TflfPHI, COMPLY, NELAST , EWALL, HW ALL, NELM, JBNUM, JTERM, 

1K,DD,CC,D,B,BB,CP, J, AREA, VME*AN,FBEL, AKBEL,BSIGN,DAMPER, SPRINGK 

COMMON G1,G2,RH0WALL,ZX,ZY 
FRICT=FPVCf J) 

AKB=AKBPVC(J) 

AKC=AKCPVCCJ) 

AKPVC=AKB-AKC 
DD(1,1, J)=FRICT 
DD(l,2,J)=(0.f D.) 

DDC2,1, J)=(0.,0.) 

DD(2,2,J)=Cl.f Q.) 

CCC1,1, J) = (Ow 0.) 

CCC2 r l, JJaAKPVC 

RETURN 

END 
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SUBROUTINE FOUR 

C THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE TRANSFER MATRIX FOR 

C A SIDE BRANCH PUL.SER 

COMPLEX ETA, ENUM , DENOM , TRANS , B , D, X I , ARG , R J , GAMMA, ZC , COSHG, SINHG, DD 
l,CC,TCnTH,TCSCH,BB,Zl,ZarXl,Xa,X3,X»f,QUADRAT,ALPHAP,BETAP,ALPHADP, 
SBETADP,GOFS,SLCW,COSHCW,SlNHCWrCOEFF, ALPHA 1, ALPHAS 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION ITYPEC1S), JTERMf 15), K( 15,15 ),EL( 15), RADIUS CIS), RBUBC 15), 
1FPVC(15), AKBPVCQ5) , AKCP VC ( IS ) , PARLE_N C 15 , 5 ) , P ARR ADX1-5 '•5-)-,-FRE-QC-aoa-)- 
- DIMENSION BTS, 1), D (?, SI, SIZE (aba ), SIZEDB (Sno), ANGLE ( 500), FBEL CIS), 
1C0MPLY CIS), AKBELC ISO, BSIGNC10), DAMPER C IS ),SPRINGK ( 15 ) 

DIMENSION DDC5,a,15) , CC C2 , 1 , IS) , BB (E , 1 , 15 ) , R ADSEC (200 ) , AREA Q5) 
COMMUN PI,RADIUS,OMEGA,ENU,ETA,EL,CO,RHOQ,THERMK,RBUB,CPCV,PO, vise 
1 , FP VC , AKBP VC , AKC PVC , NPAR, PARR AD, PARLEN , BULKMOD , RHOL IQ , NGAS, G ASMW 
COMMON GAMGAS, G, TO , PHI , COMPL Y , NELAST , EWA LL , HW ALL, NELM , JBNUM , JTERM, 
1K,DD,CC.,D,B,BB,CP, J, AREA, V MEAN, FBEL, A KBEL,BS I GN, DAMPER, SPRINGK 
COMMON G1,GE,RH0WALL,ZX,ZY 
DD(1,1,J) = U.,0.) 

00(1,2, J) = C0.,0. ) 

DD(a,l, J) = C 0 ., 0 . ) 

DDCE,a, J) = C 1 ., 0 . ) 

CC(1,1,J) — CO.,0.) 

ccca,i,j)=ci.,o.) 

RETURN 

END 
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SUBROUTINE FIVE 

THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE TRANSFER MATRICES 
FOR THE LONGITUDINAL RIGID BODY MOTION OF A SIMPLE LINE OF LENGTH 
L. THE MOTION IS A RESULT OF VELOCITY EXCITATION. 

COMPLEX ETA,ENUM,DENOM,TRANS* B,D,XI, ARG,RJ,GAMMA,ZC,COSHG,SINHG,DD 
1,CC, TCOTH,TCSCH,BB, Z1 ,Z2, XI, X2,X3, XM-, QUADRAT, ALPHAP, BET AP, ALPHADP , 
2BETADP,G0FS,SLCW,C0SHCW,SINHCN,C0EFF, ALPHA!, ALPHA5 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION ITYPE(IS), JTERM 115 ) , K ( 15 , 1 S ) , EL ( 15 J , RADIUS ( IS) , RBUB (15 ) , 
1FPVC(15) ,AKBPVC(15) , AKCPVCC15) , P ARLEN ( 1 5 , 5 ) , PARRAD ( 15 , 5 ) , FREQC200 ) 
DIMENSION B(2,l) , D(2,2) , SIZE C 200) ,SIZEDB( 2 Qn), ANGLEC200) ,FBEL(15) , 
1C0MPLYC15) , AKBELU5) r BSIGN(IO) , DAMPER C 15) ,SPRINGK(15) 

DIMENSION DDO,2,15),CC(2,l,15) , BB ( 2 , 1 , 15 ) , R ADSEC (200 ) , ARE A (15 ) 
COMMON PI , RADIUS, OMEGA, ENU, ETA, EL, CO, RHOO, THERMK, RBUB, CPC V,PO , VISC 
1,FPVC, AKBPVC, AKCPVC,NPAR,PARRAD,PARLEN,BULKMOD,RHOLIQ,NGAS,GASMW 
COMMON GAMGAS»G, TO, PHI, COMPLY, NE LAST, EWALL,HWALL,NELM,JBNUM,JTERM, 
J K,DD,CC, D,B,B6, CP, J, AREA, V MEAN, FBEL,AKBEL,8SIGN, DAMPER, SPRINGK 
COMMON G1,G2,RH0WALL,ZX,ZY 
CALL SPEED 

REYN0LD=VMEAN*2.*RADIUS(J)/(ENU*12.) 

RT = 1T'T.*(2.*ENU*0.Q055*REYNGLD**0.85)/(RADIUS(J)**2) 

Z1=COMEGA*EL(J)*ETA)/CO 

Z2=CS0RT(( 1. ,0. )+RT/ IOMEGA* ETA) ) 

XI=SQRT(OMEGA/ENU)*CSQRT(-ETA) 

ARG=XI*RADIUS(J)/12. 

CALL JIQRJD (ARG,RJ) 

GAMMA=ETA*OMEGA*ELCJJ/CC0*CSQRTC(l. ,0.)-2.*RJ/ARG))+REALCZl*Z2) 
7C=RHOO*CO/(CSQRT(( 1. , 0 . ) -2 . *R J/ ARG) ) 

COSHG= CCEXPC GAMMA )+CEXP (-GAMMA) )/2. 

S I NHG=(CEXP (GAMMA )-CEXP (“GAMMA ))/2. 

DD(1,1, J )=COSHG 

DD(1,2, J)=-ZC*SINHG/AREA(J) 

DD(?,1, J)=-AREA(J)*SINHG/ZC 

DD(2,2, J)=COSHG 

CC(1 , 1, J)=ZC*SINHG 

CC(2,1, J)=AREA(J)*((1.,U.)-C0SHG) 

RETURN 

END 
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SUBROUTINE SIX 

THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE TRANSFER MATRIX FOR 
AN ARBITRARY NUMBER OF LINES IN PARALLEL ALL OF WHICH EMANATE FROM 
AND CONVERGE ON COMMON INPUT-OUTPUT POINTS. 

COMPLEX ETA , ENIJM, DENOM, TRANS, B,D, XI ,aPG,RJ, GAMMA, ZC, COSHG,SINHG,DD 
1,CC« TCOTH,TC$CH,BBr Zl, 22, Xl,X?,X3, XY, QUADRAT/ ALPHAP, BETA P, ALPHADP, 
2BETADP,G0FS,SLCW,C0SHCW,SINHCW, COEFF/ ALPHAlr ALPHA2 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION ITYPE(lS)rJTERM.C 15 J. , K.C 1 5., 1 5-1, E L ( 1-5-)-, R A D-I-U S C IS ) /R B U SC IS) , 
1F’PVC(15) /AKBPVC(i5)/AKCPVC(i5)/PARLENflR,5),PARRAD(l5/5)/FREQ(20Q) 
DIMENSION BC2/1)/ D (2 , 2 ) / S I ZE ( 200 ) / SIZEDB (200 ) , ANGLE ( 200 ) , FBELC IS ) , 
1C0MPLYC15}/ AKBEL(15)/B3IGN(10),DAMPERC15),SPRINGKC15) 

DIMENSION DD(2,2, 15) , CC ( 2 , 1 , 1 5 ) , BB ( 2 , 1 , 15 ) , R ADSEC (200 ) , ARE A ( 1 5 ) 
COMMON PI, RADIUS, OMEGA, ENU, ETA, EL/ CO, RHOG, THERMK, RBUB, CPC V,PQ, VISC 
1 , FPVC, AKBPVC, AKCPVC, NP4R, P4RRAD, PARLEN,BULKMOD, RHOLIQ/NGAS, GASMW 
COMMON GAMGAS, G, TO , PHI, COMPLY, NELAST , EWALL, HW ALL, NELM, JBNUM, JTERM, 
1K,DD,CC,D,8,BB,CP, J, AREA, VMEAN,FBEL, AKBEL, BSIGN, 0 A MPER , SPRINGK 
COMMON G1,G2,RH0WALL/ ZX,ZY 
DO 5 1=1 , NPAR 
CALL SPEED 

CO=CO/SQRTCl.+BIJLKMOD*2.*PARRADr J, I)/f ,*EWALL*HWALD) 
SGJFT=PI*PARRADCJf I)**2/ltf . 

XI=SQRT(OMEGA/ENU)*CSORTf-ETA) 

APG=XI*PARRADCJ, I)/12. 

CALL J10RJ0(ARG/RJ) 

GAMMA=ETA*OMEGA*PARLEN(J,I)/(CD*CSQRT( Cl.,0. )-2.*Rj/ARG)) 
ZC=RHOD*CO/(CSQRT(Q. , 0 . )-2 . *R J/4RG) ) 
C0SHG=(CEXP(GAMMA)+CEXPOGAMMA))/2. 
SINHG=(CEXP(GAMMA)-CEXP(-GAMMA))/2. 
TC0TH=TC0TH+S0FT*C0SHG/(SINHG*ZC) 

5 TCSCH=TCSCH+SQFT/(ZC*SINHG) 

DU(1,1,J)=TC0TH/TCSCH 

DD(1,2, J)=-1./TCSCH 

DD(2,1, J)=TCSCH-TCOTH**?/TCSCH 

DD(2,2, J)=TCOTH/TCSCH 

RETURN 

END 
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SUBROUTINE SEVEN 

THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE TRANSFER MATRIX FOR 
A BELLOWS WITH A LUMPED COMPLIANCE REPRESENTING 1 RAPPED GAS UNDER 
THE LINER PLUS A PERIODIC VOLUMETRIC CHANGE DUE TO RELATIVE MOTION 
OF THE ENDS OF THE BELLOWS 

COMPLEX ETA,ENllM,DENOM, TRANS, B,D, XI, ARG,R J, GAMMA , ZC , COSHG, SINHG , DD 
1,CC, TCC1TH,TCSCH,BB,Z1,Z2, XI, X2,X3,X>+, QUADRAT , ALPHAP, BE TAP, ALPHADP, 
2BETADP,G0FS,SLCW,C0SHCW,SINHCW, COEFFr ALPHA1, ALPHAS 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION ITYPEQ5), JTERM(15) ,KC15,15) ,EL(15), RADIUS (15), RBUBQ5), 
lFPVC(15),AKBPVC(15),AKCPVC(iS),PARLENC15,S),PARRAD(lS,5),FREQC800) 
DIMENSION B(2,l),D(2,2),SIZE(20a),SIZEDB(2no) , ANGLEC 200 ) , EbELCl 5 ) , 
1C0MPLY( 15), AKBELC15),.BSIGN CIO), DAMPER (IS ),SPRINGK( 15) 

DIMENSION DDC2,2,lS),CCCe,l,15),BB(2,l,15),RADSECC200),AREACl5) 
COMMON PI,RADIUS,OMEGA,ENU,ETA,EL,CO,RHOO, thermk,rbub,cpcv,po, vise 
1,FPVC,AKBPVC, AKCPVC,NPAR,PARRAD,PARLEN,8ULKM0D,RH0LIQ,NGAS,GASMW 
COMMON GAMGAS,G,TO,PHI, COMPLY, NEL AST, EW ALL, HWALL,NELM, JBNUM, JTERM, 
1K,00,CC,D,B,BB,CP, J, AREA, VMEAN,F8EL, AKBEL,BSIGN, DAMPER, SPR1NGK 
COMMON G1,G2,RH0WALL,ZX,ZY 
FRICT=FBEL(J) 

C=COMPLYCJ) 

AKB=AK8EL(J) 

DDC1,1,J)=FRICT 

DD(1,2,J)=(0.,0») 

DD(2, L, J)=-C*ETA*OMEGA 
DD(2,2,J)=C1.,0.) 

CC(1,1, J)=CO. ,0.) 

CC (2 , 1 , J)=AKB 

RETURN 

END 
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SUBROUTINE EIGHT 

THIS SUBROUTINE CALCULATES THE ELEMENTS OF 1 HE TRANSFER MATRICES 
FOR THE AXIAL MOTION OF A VERTICAL LINE WITH MOUNTING STIFFNESS 
ON THE TWO ADJACENT HORIZONTAL LINES. THE TWO STIFFNESS TERMS ARE 
LUMPED INTO ONE EFFECTIVE SPRING IN PARALLEL WITH ONE EFFECTIVE 
VISCOUS DAMPER. THE FORCING FUNCTION IS AN ACCELERATION APPLIED 
TO THE SPRING-DAMPER SUPPORT PLUS A PASSIVE INERTIAL LOADING 
COMPLEX -ETA, ENUM, DENOM* TRAMS, B, D, XI ,ARG,RJ, GAMMA, ZC , COSHG , SINHG , DD 
1, CC,TC0TH,TCSCH,BB,Z1,Z2, XI, X2,X3, XI, QUADRAT, ALP HAP, BETA P„, ALP HA D„P, 
2 B L-T-A-D P; G 0 F S , S L C W, C O S H C w , STN FTC' W , C 0 E F F , A L P H A 1 , A L P ri A 2 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION I TYPE (15) , JTERM ( IS ) , K ( 1 5 , 1 5 ) , EL ( 15 ) , R ADI US C 15 ) , RBUB ( 15 ) , 
1FPVCC15), AK6PVCC15), AKCPVC(15) , P ARLEN C 1 5 , 5 ) , P ARR AD ( 1 5 , 5 ) , FREQ ( 20 D ) 
DIMENSION BC2, 1) ,D(2,2) , SIZE ( 200 ) , SI ZEDB( 200 ) , ANGLE (200 ) ,FBEL(15) , 
lCUMpLY(lS) , AKBELQS) , BSIGN(IO) , DAMPER CIS) , SPRINGK(IS) 

DIMENSION DD(2,2,1S) ,CC(2,1,15) ,8B(2,1,15) ,RADSEC(20n) , AREA(IS) 
COMMON PI,RADIIJS,OMEGA,ENU,ETA,EL,CO,RHQO,THERMK,RBUB,CPCV,PO,VISC 
1 , FP VC , A k BP VC ,AKCPVC,NPAR,PARKAD,PARLEN, BULK MOD, RHOLIQ,NGAS,GASMW 
COMMON GANG AS, G, TO, PHI, COMPLY, NELAST,EWALL,HWALL,NELM,JBNUM, JTERM, 
IK, DD,CC,D,B,BB, CP, J, AREA, VMEaN,FBEL,AKB£L,BSIGN, DAMPER, SPRTNGK 
COMMON G1,G2,KH0WALL,ZX,ZY 
CALL SPEED 

REYN0LD=VMEAN*g,*R4DlUS(J)/(ENU*12.) 

RT = 14-4- .*(2.*ENU*0. L10S5*REY MOLD* *U. 85) /(RADIUS ( J)**2) 

ZJ=(OMEGA*EL(J)*ETA)/CO 

72=CSQRT((1. ,0.)+RT/(0MEGA*ETA) ) 

Xl=SQRT(OMEGA/ENU)*CSQRT(-ETA) 

ARG=XI*RADIUS(J)/ig. 

CALL J10RJO(ARG,RJ) 

CAMMA = ETA*OMEGA*EL(J)/(Crj*CSURT( (1. , 0 . )- 2 .*RJ/ARG))+REAL( J.i *Zd') 

ZC = RH00*C0/(C3QRT( (!.,(!.) -2. *RJ/ARG)) 

C0SHG=(CEXPC GAMMA) +CEXP(-GAMMA))/2. 

S I NhG=s(CEXP( GAMMA )-C£XP (-GAMMA ))/2. 
AMASS=RHOO*CAREA(J-1)*EL(J-1)+AREA(J+1)*EL( J+O) 

QUADRAT = -AMASS*0MEGA**2 + ETA*C)MEGA*DAMPER( JJ + SPRINGK ( J) 
ALPHAP=(ETA*OHEGA*DAMPER( J)+SPR I NGK(J) )*ZC*SlNHG/( QUADRA T*FTA*OMEG 
IA) 

PETAP=AREA(J)*ErA*OMEGA*ZC*SlNHG/QUADRAT 

ALPHADP=(E'fA*OMEGA*DAMPERC J) + SPRINGK( J) )*AREA( J)*( (1. ,0. ) -COSHG) /( 
1FTA*0MEGA*QUADRAT) 

RE1ADP=(APEA(J)**2)*ETA*0MEGA*( (1. ,0. ) -COSHG) /QUADRAT 
DD(l,i, J)=CCOSHG+0FTAP)/C1.+BETAP; 

nod, 2 , j)=-zc*sinhg/(aheacj)*(i.+betap)) 

DU (2,1, J)=-AREA(J)*SINHG/ZC+BETADP*C(1. , 0 .) “COSHG) / (1 . +BETAP) 
00(2,2, J)=COSHG + BET A(;P*ZC*S1NHG/( AREA CJ)*(1. + BETAP) ) 

CCC1, l, J)=ALPHAP/(1.+BETAP) 

CC(2, 1, J)=ALPHADP/(1.+6ETAP) 

RETURN 

END 
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SUBROUTINE NIN.E 

THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE TRANSFER MATRIX 
FOR A LINE WITH FORCED CHANGES IN LINE LENGTH 

COMPLEX ETA,ENUM,DENOM,TRANS, B,D f XI , ARG, R J, GAMMA, ZC, COSHG, SI NHG, DO 
l,CC,TCOTH,TCSCH,BB,Zl,Z2,X] , X? , X^ , X4- , QUADRAT , ALPHAP , BE T AP, ALPHADP, 
2BETADP,GOFS,SLCW,COSHCW,SINHCW,COEFF, ALPHA1, ALPHAS 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION I TYPE (15), JTERM ( 15 ) , K ( 15 , 15 ) , EL ( 15 ) , R AD IUS( 15 ) , R8UB(15) , 

1FPVC(15), AKBPVC(15),AKCPVC(15),PARLEN(15,5),PARRAD(15,5),FREQ(200) 

DIMENSION B(2,1),D(2,2),SIZE(2Q0) , SIZEDBC2QG) , ANGLE ( 200) , FUEL ( 1 5) , 
1C0MPLYC15) ,AKBEL(1S) ,BSIGN(1U) , DAMPER (15) ,SPRINGK(15) 

DIMENSION DD(£,2,15),CC(2,1,15),BB(2,1,15) , R ADSEC ( 200 ) , AREA ( 15 ) 
COMMON PI, RADIUS, OMEGA, ENU, ETA, EL, CO, RHOO , THERMK, RBUB, CPCV, PO , VI SC 

1,FPVC, AKBPVC, akcpvc,npar,parrad,parlen,bulkmod,rholiq, NGAS,GASMW 

COMMON G AMG AS, G, TO, PHI , COMPLY, NEL AST, EWALL,H WALL, NELM, JBNUM,JTERM, 
1K,DD,CC,D,B,BB,CP, J, AREA, VME AN, FBEL , AKBEL , BSIGN , DAMPER, SPR1NGK 
COMMON G1,G2,RH0WALL,ZX,ZY , 

G0FS=G1+G2*ETA 

CW=SQRT((EWALL*G)/(RH0WALL*12. )) 

CALL SPEED 

REYN0LD=VMEAN*2.*RADIUS(J)/(ENU*I2.) 

RTslf^.*(2.*ENU*0.005S*REYN0LD**0.85)/(RADIUS( J)**2) 

Z1=(0MEGA*EL(J)*ETA)/C0 

Z2=CS0RT( (1. , 0. )+RT/(OMEGA*ETA) ) 

XI=SQRT(OMEGA/ENU)*CSQRT(-ETA) 

ARG=XI*RA0IUS(J)/12. 

CALL J10RJQ(ARG,RJ) 

GAMMA=ETA*OMEGA*EL( J)/(CO*CSQRT ( ( I . , 0 . ) -2 . *R J/ ARG) ) +Rb AL ( 7. 1*Z2 ) 
ZC=RH00*C0/(CSQRT((1. ,0.)-2.*RJ/ARG) ) 

COSHG=(CEXP( GAMMA) +CEXP (-GAMMA) )/2. 

SINHG=(CEXP (GAMMA) -CEXP (-GAMMA) )/2. 

DD(1,1, J)=COSHG 

DD(1,2, J)=-ZC*SINHG/AREA(J) 

DU(2,I,J)=-AREA(J)*SINHG/ZC 

DD(2,2, J)=COSHG 

SLCW=ETA*OMEGA*EL(J)/CW 

C0SHCw=(CtXP(SLCW)+CEXP(-SLCW))/2. 

SINHCWs ( CEXP (SLCW) -CEXP (-SLCW) )/2. 

C0EFF=(-UMEGA**2/C0**2- ( GAMMA/EL ( J) )**2)/ (-0MEGA**2/CW**2- (GaMMA/E 

1L(J))**5) 

ALPHA1=(RH00*C0**2/(ETA*0MEGA) )*( (eta*omega*sinhcw/cw-gahma*sinhg/ 

1EL(J))+(GOFS-COSHCW)*(COSHCW"COSHG)*ETA*OMEGA/(CW*SINHCW)) 

ALPHA2 = - ( ( COSHC W-COSHG) + ( GOFS-COSHC W) * CS I NHCW-ET A*OMEGA*EL( J)*SINH 
1G/(CW*GAMMA) )/SINHCW) 

CC(1,1, J)=C0EFF*ALPHA1 
CC(2,.l , J) = C0EFF*ALPHA2*AREA( J) 

RETURN 

END 
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SUBROUTINE TEN 

THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE TRANSFER MATRIX 
FOR THE AXIAL MOTION OF A VERTICAL LINE WITH MOUNTING STIFFNESS. 
THE STIFFNESS IS REPRESENTED BY THE IMPEDANCE PRESENTED TO THE 
LINE BY THE VFhICLE STRUCTURE. 

THIS SUBROUTINE CAN BE USED ONLY WHEN ANOTHER SOURCE OF EXCITATION 
IS PRESENT IN THE FEEOLINE , I.E. PULSER OR VERTICAL LINE VELOCITY 
EXCITATION. THE FUNCTIONAL FORM OF THE STRUCTURAL IMPEDANCE IS 
THAT OF AN EQUIVALENT SPRING AND DASHPOT IN PARALLEL. 

ZX AND Z-Y -ARE- THE -RE-AL -AND- IM AG-PNA-RY' PARTS OF ZSTRUCT.' 

COMPLEX ETA,ENUM,DENOM, TRANS, B, D , XI , ARG, RJ, GAMMA, ZC, COSHG, SINHG, DD 
If CC, TCOTH, TCSCH,Bb,'Zl, Z8, XI , X 5, X3,X>f , QUADRAT, ALPH AP, BETA P,ALPHADP, 
2BETADP , GOFS, SLCW , COS HCW, SI NHCW , COEFF, ALP HA1 , ALPHAS 
CUMPLEX ZSTRUC T r ALPHA, BETA 

DIMENSION ITYPEC15), JTERM C 15) , K( IS , 15 ) , EL (15) , RADIUS (15 ) , RBUB ( 15 ) , 
lFPVC(lS) ,AKBPVC(15), AKCPVCf IS) , PARLEN (i 5 , 5) , PARRADC15, S) , FREQ(HOO) 
DIMENSION B(3, 1), DC 8, 8), S I ZE ( 8 (3 0 ) , S I ZEDB ( 800 ) , ANGLE ( 8 Q 0 ) , FBEL ( 15 ) , 
1C0MPLY(15),AKBEL(15),B3IGN (10), DAMPER (15), SPRINGKC15) 

DIMENSION DD(8,8,15) ,CC(8, 1 , 1 5 ) , BR ( 8 , 1 , IS ) , R ADSEC ( 8 OO ) , AREA (15) 
COMMON PI, RADIUS, OMEGA, ENU, ETA, EL, CO, RHO[),THERMK,RBUR,CPCV,PO,V ISC 
1 ,FPVC, AKBPVC, AKCPVC, NPAR,PARRAD, PARLEN, BULKMOD,RHOL IQ, NGAS,GASMW 
COMMON GAMGAS, G, TO, PHI, COMPLY ,NEL AST, E WALL, H WALL, NELM,JBNUM,J TERM, 
1K,DO,CC,0,B,BB,CP, J, AREA, VMt AN, FBEL,AKBEL,BS IGN, DAMPER, SPR I NGK 
COMMON G1,GS,KHQWALL,ZX,ZY 
CALL SPEED 

PEYNOL0=VMEAN*a.*RADIUS( J)/(ENU*18. ) . 

RT = l l f i l.*(a.*ENU*li.0055*REYN0LD**rj.B5)/(RADIUS(J)**5) 

71=(0MEGA*EL( J)*ETA)/CO 
7d=CSQRT( (1. ,0. )+RT/(OMEGA*ETA) ) 

XI=SQRT(OMFGA/ENlJ)*CSQRT(-ETA) 

ARG=X1*RADIUS(J)/1?. 

CALL J10RJ0(ARG,RJ) 

GAMMA=ETA*OMEGA*EL( J)/(C0*CSQRT((l.,n.)-8.*RJ/ARG))+REALCZl*Z8) 
ZC=RHOO*CO/(CSQRT((1. ,0. ) -8 . *R J/ A RG > ) 
C0SHG=(CEXP(GAMMA)+CEXP(«GAMMA))/2. 
SIMHG=(CEXP(GAMMA)-CEXP(-GAMMA))/2. 

AMASS=RHOn*( AREA ( J-l ) *EL C J- 1) + ARE A ( J + l ) *EL C J+l ) ) 

ZSTRuCT=ZX-ET A*ZY/OMEGA 

ALPHArAREA( J)*ZC*SINHG/ (ZSTRUC T+AMASS*ETA*OMEGA) 

BET As (l.-COSHG )*( AREA (J)**8)/ (ZSTRUC T +AMASS*ETA*OMEGA ) 

00(1, 1, J)=(C0SHG+ALPHA)/(1.+ALPHAJ 
DDll, 2, J)=-ZC*SINHG/(AKEA(J)*(1.+ALPHA)) 

DO (8,1 , J)=-A«EA ( J)*SINHG/ZC + nETA*(l.-COSHG)/(,l. +ALPHA) 

DD( 8,8, J)=C0SHG+ZC*SINHG*BETA/(AREA(J)*(1.+ALPHA) ) 

RETURN 

FND 
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SUBROUTINE ELEVEN 

THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE TRANSFER 
MATRIX FOR A COMPLEX SIOE BRANCH. THIS MODEL ASSUMES A 
LAMINAR FLOH IN THE SIDE BRANCH* 

COMPLFX ETA,ENUM, DENOM, TRANS, 8,D, XI ,ARG,RJ, GAMMA, ZC, COSHG, SINHG, DD 
1, CC,TCOTH,TCSCH,BB,Zl,Z2, Xl ,X2,X3,X4, QUADRAT, ALPHAP, BETAP, ALPHADP, 
2BETADP,G0FS,SLCW,C0SHCW,SINHCW,C0EFF, ALPHA1,ALPHA8 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION I TYPE (15), JTERM(IS)., K(15,15) , ELY 15 ) , R AD IUS ( 1 5 ) , RBUB(IS), 
lFPVCa5),AKBPVC(15),AKCPVCC15),PARLENQS,5),PARRADCl5,5),FREQ(200) 
DIMENSION B(2,1),DC2,2:) ,SIZEC2Q0 ),SIZEDB (200), ANGLE (200), FBELC15), 
1C0MPLYQS),AKBEL(1S),BSIGN(10) , DAMPER (15) , SPRINGKQS) 

DIMENSION DD(2,2,15),CC(2,l,15),BB(2,l,15),RADSEC(20fl),AREA(15) 
COMMON PI,RADIUS,OMEGA,ENU,ETA,EL,CO,RHOO, THERMK,RBUB,CPCV,PO, VISC 

1,FPVC, AKBPVC, AKCPVC,NPAR,PARRAD,PARLEN,BULKMQO,RHOLlG,NGAS,GASMw 

COMMON GAMGAS,G, TO, PHI, COMPLY, NELAST,EWALL,HWALL,NELM,JBNUM,JTERM, 
1K,DD,CC,D,B,BB,CP, J, AREA, VME AN, FBEL,AKBEL, BS I GN, DAMPER, SPRINGK 
COMMON G1,G2,RH0WALL,ZX,ZY 
COMMON BRC0MP,BRL,BRDIAM 
COMPLEX S 
S=ETA*OMEGA 
10 CALL SPEED 

AREAB=( PI/*. )*(BRDI AM/12. )**2 
BRI=RH00*BRL/AREA6 

BRR=128.*VISC*BRL/(PI*(BRDIAM/12. )***.) 

00(1,1, J)=Cl.,P. ) 

DD(1,2, J)=(0. ,0.) 

DD(2,1, J)=-S*BRC0MP/(8RI*BRC0MP*S**2+BRR*BRC0MP*S +1. ) 
DD(2,2,J)=(1.,G.) 

RETURN 

END 
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SUBROUTINE SPEED 

THIS SUBROUTINE CALCULATES THE SPEED OF SOUND FOR THE FLUID IN THE 
FEEDLINE. THIS FLUID MAY BE A LIQUID OR A HOMOGENEOUS MIXTURE OF 
A LIQUID AND AN ULLAGE GAS 

COMPLEX ETA,ENUM,DENOM, TRANS, B,D, XI, ARG,RJ, GAMMA, ZC , COSHG, SINHG, DD 
1,CC,TC0TH,TCSCH, BB,Z1,Z2, Xl , X2,X3, X4, QUADRAT, ALPHAP, BETAP , ALPHADP, 
2BETADP,G0FS,SLCW,C0SHCW,SINHCW.,C0EFF.,.ALPHA1,-ALPHA2- 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION I TYPE (15), JTERM C 15 ) , K C 15 , 15) , EL (15 ) , RADIUS(IS) r RBUBC15) , 
1FPVC(15) , AKBPVCC15) ,AKCPVC(15),PARLEN(15,5), PARRAD (15 , 5) r FREQC200) 
DIMENSION 6(2,1),D(2,2),SIZE(200),SIZEDB(2QO),ANGLE(200),FBEL(15), 
1C0MPLYC15),AKBEL(1S),BSIGN(10),DAMPER(15),SPRINGK(15) 

DIMENSION DD(2,2,15),CC(2,1,15),BB(2,1,1S),RADSECC200),AREA(1S) 
COMMON PI, RADIUS, OMEGA, ENU, ETA, EL, CO, RHOO,THERMK,RBUB,CPCV,PO, VISC 
1,FPVC,AKBPVC,AKCPVC, NPAR, P ARRAD, PARLEN, BULKMOD, RHOLIQ, NGAS, GASMW 
COMMON GAMGAS, G, TO , PHI , COMPLY, NELAST, EW ALL, HWALL, NELM,JBNUM, JTERM, 
1K,DD,CC,D,B,BB,CP,J,AREA, VMEAN, F8EL, AKBEL, BSIGN, DAMPER, SPRINGK 
COMMON G1,G2,RH0WALL,ZX,ZY 
CLIQUID=SQRT (BULKMOD /RHOLIQ) 

IF(NGAS.EQ.1)10,20 

10 RH0GAS=144.*P0*GASMW/(l545,*(TD + H'bD.)*G) 
CGAS=SQRT(GAMGAS*G*1545.*(T0+4bO.)/GASMW) 

El =RHOGAS+PHI* RHOLIQ 

E2=RHOGAS/(RHOLIQ*CLIQUID**2*E1)+(PHI*RHOLIQ)/(E1*RHOGAS*CGAS**2) 

E3-(1.+PHI)*RHOGAS*RHOLIQ 

C0=SQRT(E1/(E2*£3)) 

RH00=E3/El 
ENU = VI SC/RHOO 
GO TO 30 
20 C0=CLIQUID 
RHOO=RHOLIQ 
tNU=V ISC/RHOO 
30 IF(NELAST.EQ.l) l f0,50 

40 C0=C0/SQRT(1.+C0**2*RH00*2.*(RADIUS(J)+HWALL)/(144.*EWALL*HWALL)) 
50 RETURN 
END 
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' SUBROUTINE DM AT 

C THIS SUBROUTINE CALCULATES MATRIX D IN THE EXPRESSION P=DQ+VB 

COMPLEX £TA, ENUM>QENOM> TRANS, 8, D, XI , ARG,.RJ.,. GAMMA, ZC , COSHG, SINHG, DD 
1, CC, TCOTH, TCSCH, BB, Z 1, Z5, XI, X P, X3, Xf, QUADRAT, ALPHAP, BETA P,ALPHADP, 
2BETADP, GOES, SLCW,COSHCW,SINHCN,COEFF, ALPHAl, ALPHA 2- 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION ITYPE(IS) , JTERM(iS).,K(15,i5) ,EL.(15) ,RADIUS(15),RBUB(15) , 
1FPVCC15), AKBPVCC15), AKCPVC ( 15 ) , P ARLEN ( 15 ,3 >,-P.ARRAD(15.,5) , FREQ (200) 
DIMENSION 8(2,1) , D ( 2 , 2 ) , SI ZE C 2 00 ) , SI ZEDB ( 200 ) , ANGLE ( 200 ) , FBEL( 15 ) , 
1C0MPLYC15), AKBELC15) , BSIGNC 10) , DAMPER (15) , SPRINGKC15) 

DIMENSION DD(2,2,15) ,CC('2,1 , 15) , BB ( 2 , 1 , 15 ) , R ADSEC ( 2.00 ) , AREA (15) 
COMMON PI, RADIUS, OMEGA, ENU, ETA, EL, CO, RHOO,THERMK,RBU.B,CPCV,PO, VISC 
1,FPVC, akbpvc,.akcpvc,npar,parrad.,parlen,bulkmod,rholiq,ngas,gasmw 
COMMON GAMGAS,G,T.O,PHI,COMPLY,NELAST,EvJALL,HWALL,NELM,.JBNUM, jterm, 
1K,DD,CC, D, B,BB, CP, J, AREA, VME AN, FBEL,AKBEL,BS IGN, DAMPER, SPRINGK 
COMMON G1,G2,RH0WALL,ZX,ZY 
DO 10 M=l,2 
DO 10 N=l,2. 

10 D(M,N)=DD(M,N,1) 

IF(NELM-1)S0,20,30 
20 RETURN 
30 DO + 0 J=2 , NELM 

X1=DD(1,1, J)*D(1,1)+DD(1,2, J)*D(2,1) 

X2=DD(1,1, J)*D( 1,2) +00(1,2, J)*D(2,2) 

X3=DD(2,1, J)*D (1,1) +00(2,2, J) *0(2,1) 

X+=D0C2,1, J)*D(1,2)+DD(2,2, J)*D(2,2) 

0(1,1)=X1 
D(1,2)=X2 
0(2, 1)=X3 
+ 0 DC2,2) = X + 

RETURN 
50 PkINT bO 
CALL EXIT 

bO FORMAT (28HOERROR IN NELM SPECIFICATION) 

END 
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SUBROUTINE BMAT 

C THIS SUBROUTINE CALCULATES MATRIX B IN THE EXPRESSION P=DQ+VB 

COMPLEX ETA, ENUM,DENOM, TRANS, B,D, XI , ARG,ftJ, GAMMA, ZC , COSHG , S INHG, DD 
l,CC,TC0TH,TCSCH,BB,Zl,Z2,Xl,X2,X3,X‘f, QUAD RAT, ALPHAP,BETAP,ALPHADP, 
26ETADP,G0FS,SLCW,C0SHCW,SINHCW,C0EFF, ALPHA1, ALPHAS 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION ITYPE(IS) , JTERM(is) ,K(1S,15) ,EL(15) , RADIUS (IS) , RBUB CIS) , 
1FPVC(15),AKBPVC(15),AKCPVC(15),PARLEN(15,5), PARR AD. (1.5. , 5.).,. F.R E Q.(-2-0 0 ) 
DIMENSION B( 2,1), DC 2, 2), SIZE (200) , SIZEDB (200 ), ANGLE (200 ), FBELC15 ) , 
1CGMPLYUS), AKBELQ5) ,BS IGN (10), DAMPER (15) ,SPRINGK(15) 

DIMENSION DD(2,2,15),CC(2,1,15),BB(2,1,15),RADSEC(2G0),AREA(15) 
COMMON P 1, RADIUS, OMEGA, ENU, ETA, EL, CO, RH0O,THERMK,RBUB, CPC V,PO, VISC 

1,FPVC, akbpvc,akcpvc,npar,parrad,parlen,bulkmod,rholiq,ngas,gasmw 

COMMON GAMGA S,G, TO, PHI, COMPLY, NEL AST, Ew ALL, HWALL,NELM,J6NUM,JTERM, 
lK,Di),CC,D,B,BB,CP, J, AREA, V MEAN, FBEL,AKBEL,BS IGN, DAMPER, SPRINGK 
COMMON G1,G2,RH0WALL,ZX,ZY 
00 10 U I = 1,JBNUM 
KOUNT=JTERMCI) 

L=K ( I , KOUNT) 

DO 100 J=l, KOUNT 
IF(J-1)S,5,15 
5 DO 10 Msl, 2 
10 8B(M,1,I)=CC(M,1,L) 

GO TO 100 

15 L“K(I,K0UNT-J+1) 

X1=DDC1,1,L)*BBC1,1 , I)+D0(1,2,L)*BB(2,1, I) 

X2 = DDC2,l,L)*BB(l,i , I ) + DDC 2,2,0*88(2,1, I) 

BB(1,1,I)=X1 
BBC2, 1, I ) = X2 


100 

CONTINUE 





8(1, n 

=BSIGN(i)* 

BB( 

1,1,1) 


8(2,1) 

=BSIGN(1)* 

BB( 

2,1,1) 

1 1 n 

DU 2nn 

Msi, 

2 




DO 200 

1 = 2, 

JBNUM 



2no 

R ( M , l ) 

sBCM, 

1) + BS 

IGN 

CI)*BH(M,1,I) 

210 

RETURN 






END 
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SUBROUTINE JlORJO(Z,RJ) 

CALCULATION OF JO(Z) AND JlCZ) Z-COMPLFX 

COHPLEY Z,v!j ,J0,TERM0,lERMl,Zl,Z2,P0rQ0,PT,Ql#PH0,PHl,FZl,FZ2,RJ 
X=REAL.(Z) 

Y = AIMAGCZ) 

R=CA6S(Z) 

IF(R-18. )10U, 100,110 

100 TERMl=Z/2. 

Jl' = Z/2. 

J0=(1. ,0.) 

TERMO= ( 1 . , 0 . ) 

A = l. 

AM=15.+R 

101 TERM0=TERM0*C-CZ/2.)**2)/A**2 . 

JO=JO+TERMO 

TERMl = TERMl*(-CZ/2. )**2)/( CA + 1.)*A-)' 

Jl= Jl+TERMl 
A=A+1. 

IF(A-AM)101,101,115 

110 IF(X)111,112,112 

111 Zl=-Z 

GO TO 113 

112 71=Z 

113 PI=3. 141552b 
Z2=8.*Z1 

IF (CABS (Z2) -5000. )120, 130, 121 

120 P0=(1. ,0.)-4.S/Z2**2+3b?5./(8.*Z2**4) 

Q0="l./Z2+37. 5/72**3-55535. /(8. *72**5) 

Pl=(l. ,n.)+?.5/Z2**2-H725./(8.*Z2**f) 
Ql=3./Z2-52.5/Z2**3+727b5./(8.*Z2**5) 

GO TO 122 

121 P0=(l.,0.)-4.5/Z2**2 
Q0=-l./Z2 

Pl=(l.,0.)+?.5/Z2**2 
01=3. /Z2 

122 PH0=Zl-Rl/4. 

Phl=Zi-. 75*PI 
FZ1=2./(PI*71) 

FZ2=CS0RT(FZ1) 

AZ1=A IMAG(Zl) 

IF(ABS(AZl)-50.)llb,llb,117 
lib J0=FZ2*(P0*CC0S(PH0)-Q0*CSIN(PH0)) 

J1=FZ2*CP1*CC0S(PN1)-Q1*CSIN(PH1)) 

IF(X)114, 115,115 
111 J1=-J1 
115 RJ=J1/J0 
GO TO 115 
11? SIN=Y/ ABS(Y) 

RJ=(SIN*( 0 . ,1. )*Pl + 01)/(Pn-SIN*(0. ,1. )*Q0) 

115 RETURN 
END 
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CARD 3 


CARD 4 


ITYPE(J), J=l, A sequence of dimensionless 

NELM dimensionless numbers that describes the 

type of elements in the line 
in the order in which they 
appear beginning "with the 
element' attached to the fuel 
tank outlet. See the first 
page of program listing for 
an expanded definition. 


JTERM (J), J = 1, The number of submatrices 

JBNUM dimensionless in B (J). For example, in 

Equation (125), JBNUM =3, 
JTERM (1) = 2, JTERM (2) = 3 
and JTERM (3) = 4. 


CARD-5 K(J, M), dimensionless 

6 


4 + BJNUM 


An array that fixes the order 
of matrix multiplication in 
B(J). For example, in Equa- 
tion (125), B (3) = D s D 4 D 3 C 2 . 
Then K(3,l) = 5, K(3,2)=4, 
K(3,3) = 3 and K(3,4) = 2. 


CARD BSIGN(J), dimensionless Floating point array that de- 


(5 + JBNUM) 

CARD HERTZ I, Hz 

(6 + JBNUM) 

DELHZ1, Hz 

F1LIM, Hz 

DELHZ2, Hz 
- HERTZ F, Hz 


fines the algebraic sign pre- 
ceding each B(J). Again in 
Equation (125), BSIGN(l) = 

+ 1.0, BSIGN(2) = -1.0 and 
BSIGN(3) = -1.0. 

The initial frequency where 
response calculations begin. 

Frequency step size for fre- 
quencies between HERTZI 
and F1LIM. 

Frequency break-point which 
defines the boundary between 
regions of different grid size. 

Frequency step size for fre- 
quencies beyond F1LIM. 

The final frequency to be 
evaluated. 
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CARD SIGN, dimensionless 

<7 + JBNUM) 

TERMZ, lb -sec/ ft B 


RHOL.IQ, lb - s e c 2 / ft 4 


TO, °F 

CARD EWALL, lb /inf 

(8 + JBNUM) 

H¥ALL, in. 

Z INPUT, lb-sec/ft s 

CPCV, dimensionless 

BULKMOD, lb /ft 3 

PHI, dimensionless 

CARD VISC, lb- sec/ft 3 

(9 + JBNUM) 


Algebraic sign preceding 
the parameter K in Equa- 
tion (121). 

Terminal impedance of the 
line which reflects the com- 
bined effect of the turbopump- 
injector- engine combination. 
In its present form, TERMZ 
is a scalar resistance. 

Liquid propellant density. 

Thermal conductivity of the 
gas inside large bubbles. 

Mean propellant pressure. 

Mean propellant temperature. 

Wall modulus of elasticity. 

Wall thickness. 

Scalar flow impedance at the 
propellant tank outlet. 

The ratio of specific heats 
for the gas in the cavitation 
bubble at the propellant 
temperature. 

Bulk modulus of the pro- 
pellant. 

Vapor-liquid mass ratio for 
the bulk propellant flow. 

Propellant dynamic viscosity. 

Molecular weight of dissolved 
vapor or ullage gas for use in 
speed of sound calculation. 

Ratio of specific heats for the 
gas whose molecular weight 
is GASMW. 


GASMW, dimensionless 


GAMGAS, dimensionless 


THERMK, Btu/(hr-ft-°R) 
PO, psia 
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CP, BtuAb-°R 


YMEAN’ ft/ sec 


Specific heat at constant 
pressure for the gas in the 
cavitation bubble. 

Mean propellant feed velocity. 


The sequence and number of the remaining input cards cannot 
be stated a priori since they depend on the component structure of the 
feedline. However, the input for each line component appears se- 
quentially in the same order as the components in the line beginning 
with the component attached to the fuel outlet. Listed below are the 
inputs for each class of component. 

Simple line: EL, ft 

RADIUS, in. 


Line length. 

Line radius as measured to 
the inner surface. 


Cavitation bubble: RBUB, in. 


Bubble radius. 


RADIUS, in. Internal radius of the line 

surrounding the bubble. 


Pres sure -volume 
compensator: 


FPVC, Friction factor expressed as 

dimensionless the fraction of the inlet pres- 
sure remaining at the outlet. 


AKBPVC, ft 2 Volume change constants 

AKCPVC, ft 2 for compensator. 


Side branch pulser: 


No input required. 


Line with velocity EL, ft 

excitation: RADIUS, in. 


Same as simple line. 
Same as simple line. 


Parallel lines: NPAR, 

dimensionless 
PARLEN, ft 
PARRAD, in. 

Bellows with relative FBEL, 
motion: dimensionless 

COMPLY, ft 5 / lb 
AKBEL, ft? 


The number of parallel 
lines. 

Line length. 

Line radius. 

Bellows' friction factor, same 
definition as FPVC 
Compliance of bellows. 

Volume change constant. If there 
is no relative motion, AKBEL =0. 
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Lrine with mounting 

EL, ft 

Line length. 

stiffness (discrete 

parameter): 

RADIUS, in. 

Line radius. 


DAMPER, 

Coefficient of viscous 


lb-sec/ft 

damping. 


SPRINGK,lb /ft 

Spring constant. 

Forced changes in 

EL, ft 

Line length. 

line length: 


RADIUS, in. 

Line radius. 


Gl, G2 

Real and imaginary parts 
of G(s). 


RHOWALL, 
lb /in 3 

Wall density. 

Line with mounting 

EL, ft 

Line length. 

stiffness (impedance 

technique): 

RADIUS, in. 

Line radius. 


ZX, ZY 

Components of structural 
impedance. 

Complex side branch: 

BRL, ft 

Branch line length. 


BRDIAM, in. 

Branch radius. 


BRCOMP, ft 5 / lb 

Compliance of the complex 
side branch. 


D. 4. Deck Structure for a Normal Compilation and Multiple 

Execution on the CPC 6400 Digital Computer 

Deck structure for multiple production runs i-s shown on the 
following page. 

The anticipated run time, core memory and output linage 
information indicated on the first two control cards are for illustra- 
tion purposes only. Typically, the combined central processor/ 
print processor time for compilation, execution and listing of one 
run is approximately 25 seconds. In the case of production runs, 
time-line charges can be improved by inserting a NOLIST card in 
front of the PROGRAM CONTROL card which will delete the 
source program listing from the output. 
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H. 5. Computer Symbol List 

In order to facilitate program changes by the user, a corre- 
spondence list of the symbols used in the program is given below. 
The quantities on the left are the FORTRAN alphanumeric symbols 
in the program; on the right is the corresponding item from the 
analysis. Dummy variables, which are frequently used to facili- 
tate programming of long equations, have been omitted. 


TERMZ 

z t 

EE 


L 

RHOLIQ 

PX 

RBUB 

- _ 

R 0 

(Mass density of propellant; 
liquid phase) 

PARLEN 

— 

L p 

THERMK 

k 

PARRAD 

-- 

r P 

PO 

Po 

FBEL 


F(PV x 2 ,G) BEIj 

TO 

T 0 

COMPLY 

-- 

c 

EWALL 

Et 

AKBEL 

-- 

K B, BEL 

HWALL 

h 

FPVC 


F(PV X 2 , G) pvc 

ZINPUT 

Zi 

AKBPVC 

-- 

k b, pvc 

CPCV 

C p /c v 

ACKPVC 

-- 

K C, PVC 

BULKMOD -- 

K 

DAMPER 

-- 

b 

PHI 

$ 

SPRINGE 

-- 

k 

vise 

n 

OMEGA 


U) 

GASMW 

MW 

ENU 


V 

GAMGAS 

Y 

REYNOLD 


Nr 

CP 

C P 

RT 


Rt 

VMEAN 

v 0 

ETA 

-- 

(+i) 

RADIUS 

r 0 

XI 



5 
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JIORJO 

-- 

h/Jo 

GAMMA 


T 

ZC 

-- 

Zc 

COSHG 


cosh r 

SINHG 


. sinh r 

AREA 

-- 

A 

AMASS 

-- 

mj + m 3 

ALPHAP 

-- 

a' 

BETAP 

-- 

3 ' 

ALPHADP 


a 7 ' 

BETADP 


f 

CW 


C„ 

GOFS 


G(s) 

RHOWALL 

-- 

Pt 

SLCW 

-- 

sli/ 

COSHCW 


cosh ( s JL/ c w ) 

SINHCW 

-- 

sinh ( sLr/ c w ) 

ALPHA 1 

-- 

a i 

ALPHA2 

-- 

a 3 

CLIQUID . 




RHOGAS 


P s 

CGAS 


c s 

RHOO, 


■Po 

ENU 

-- 

V 

PRO 

-- 

$ a R 0 

D THERM 

‘ 

Di 

B THERM 


b th 

BVIS 

-- 

b VISC 

BRAD 

-- 

b RAD 

POLY 


R 

VBUB 

-- 

V BUB 

SPRING 


VPo/ v BUB 

ZSTRUCT 


Zg — Z x -iZy/ U) 

AREAB 


A b 

BRI 

-- 

I 

BRCOMP 


C 

BRR 


R 

BRL 

-- 

L b 

BRDIAM 

« — 

d b 


D-30 



D. 6. Example Problems 


This section contains five example problems which will ac- 
quaint the user with the mechanics of program setup and execution. 
The development of each problem follows the same general format: 
statement of the physical problem, generation of the system matrix 
equation, presentation of critical input data and a computer listing 
and graphical representation of the configuration response. In all 
cases, only a portion of the output is shown since the remaining 
calculations follow the same format through termination of execution. 

Problem No. 1: 


The first feedline example, shown in Figure D-l, consists of 
a rigid length of line, a sidebranch pulser and a cavitation bubble in 
series with the fuel tank and a scalar terminal impedance, R. The 
perturbation pressure at Station 1 is assumed to be negligible. This 
assumption is for illustration purposes only, since the computer pro- 
gram accepts an arbitrary scalar input impedance at Station 1. The 
physical distance between Stations 2 and 4 is assumed to be small 
compared to the total line length, L. It is required that the perturba- 
tion pressure at Station 4 be determined for oscillatory pulser in- 
puts, Qd . The first step is to set up the station numbers, component 
numbers and component matrix equations in the Laplace domain as 
indicated in the illustration. Note that the subscripts on pressure 
and flow perturbations correspond to a particular station number, 
whereas the subscripts on matrices and Ch refer to a specific 

component number. Beginning with the matrix equation for compo- 
nent number 3, successive substitution of the remaining expressions 
yields the overall line equation: 


" p 4 “ 


• 0 - 


- D3 p 2 Pl 


_p 4 /r. 


_Qi_ 


+ Q d (s) (+D 3 C 3 ) 


(D.l) 


The fundamental operations involved in the derivation of this system 
equation are common to the setup of all problems, and the user is 
required to perform these steps in order to execute the program. 

At this point, the user would make the following general ob- 
servations regarding the program input; 

(1) NELM = 3 The number of components in the line. 

(2) NEXCITE = 4 Element number providing the excitation. 
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q 4 = p 4 /r 

Mean Flow Velocity = 50 ft. /sec. 

Bubble Radius = 1.2 in. 

Line Radius =4.0 in. 


Figure D-l. Line Model For Example Problem No. 1 



(2) B = D 3 C: 3 Therefore, JBNUM = 1 

JTERM(l) = 2 {two sub- 
matrices in = B ) 
K(l,l) =3 ~ 

K(l,2) = 2 

(3) ITYPE (1) = X ITYPE(J) indicates the type of compo- 

nent corresponding to component 
ITYPE (2) = 4 number J. Refer to comment cards on 

first page program listing, 

ITYPE (3) = 2 

(4) SIGN = +1.0 The algebraic sign preceding the exci- 

tation, Q d , in Equation (D.l). 

(5) BSIGN (J) =1.0, J= JBUM=1. BSIGN (J) is the algebraic 

sign preceding the Jth term in B,= . . . The 

usefulness of this input will be demonstrated later in a 
more complex example. 

In this example, we will assume that the propellant does not contain 
any dissolved gases, and the calculated speed of sound in the pro- 
pellant will reflect wall elasticity effects. 

Using typical values for the physical parameters associated 
with the line and the propellant, a complete input package would take 
the form shown in the following listing. At this point, the reader 
should thoroughly review the setup instructions listed at the begin- 
ning of the PROGRAM CONTROL. These instructions pertain to 
program inputs between statement numbers 10 and 15, and they indi- 
cate the flow of calculations within the program. Note that, when a 
pulser is involved in the feedline, no input is required for its des- 
cription (see statement number 35). 

T-he computer inputs and outputs for this problem are shown 
on the following page. The column labeled TRANS is the magnitude 
of P 4 /Qd. Due to the flexibility of the computer code, the transfer 
function is not normalized. The reader may wish to normalize the 
output, and, in this event, a suitable normalization factor would be 
the characteristic line impedance, Z c = p 0 c 0 . The response magni-r 
tude for this configuration is presented in Figure D-2. 
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INPUT 


EXAMPLE PROBLEM NO. 1 

3 10 1 + 

1 + 2 

3 

3 2 

• + .Tn0Trt)0+n-i 

+.100000+01 +.500000+00 +.500000+02 +.50000n+00 +.b00000+02 
+ . 1110000+01 + . + 5inofJ + os + . eaounn+oi +. + bnnun-ng + .3'*?nno+02 ?ohooo+C' 3 
+ . 3 n 0 n Q n + Cl 8 +.bb0000-01 +.000000+00 +.1 + 0000 + 01 + . 1 0 8000+08 +.00000U + UU 
+ . + u8iU'0-t'5 +.OOOOUU+UO +. 000000 + 00 +.22+OUO+UO +.500000+02 
+ . 3 U 0 0 u fi + 0 2 +.+00U0U+01 
+ . ieomiti + ul +.+G0UUO + O1 

OUTPUT 


NON-DIMENSIONAL PHASE 


OMEGA-RAD/SEC 

FREQ-HZ 

TRANS-P/K 

TRANS-DB 

ANGLE-DEG 

fe.3 

1.0 

1187.728 

-23.83 

-82.. 8 

R.+ 

1.5 

1808. b3b 

-20. +0 

-83.2 

12. b 

2.0 

2+38. 83b 

-17.81 

-83.8 

15.7 

2.5 

3081.8+8 

-IS. 75 

-8 + . 5 

18.8 

3.0 

3??b. 372 

-1+.01 

-85.2 

22.0 

3.5 

+500.573 

-12. +8 

-85.1 

25.1 

+ .0 

52?+. +35 

-11.11 

-87.0 

28.3 

+ .5 

bl08.?28 

-8.83 

-88.0 

31. * 

S.O 

7020.705 

-8.53 

-88.2 

3 + . b 

5.5 

802+ • 805 

-7. + ? 

-100. + 

37.7 

b.O 

81++.22+ 

-5.33 

-101.8 

+ 0.8 

b.S 

10+0b.3b+ 

-5.21 

-103.5 

+ 4.0 

7,0 

118+b.7?3 

- + .08 

-105.3 

+ 7.1 

7.5 

13511.200 

-2.8 + 

-107.5 

50.3 

8.0 

15+58.80+ 

-1.7? 

-110.1 

53.+ 

8.5 

17?bb.ll8 

-,5b 

-113.2 

5b. 5 

8.0 

20528.217 

.58 

-117.0 

55.7 

8.5 

23855.38b 

2.00 

-121.8 

b2.8 

10.0 

278+7.5+3 

3.3 + 

-128.0 

bfa.G 

10.5 

3251S.+b8 

+ .58 

-135.8 

58.1 

11.0 

37580.835 

5.85 

-1 + 5.1 

?2.3 

11.5 

+2238. 5b3 

5.8b 

-158.7 

75.4 

12.0 

+50+5.750 

7.52 

-173.2 

78.5 

12.5 

++818.255 

7. + 8 

-188.2 

81.7 

13.0 

+2185.753 

5.85 

-201. b 

8+ . 8 

13.5 

38205. +05 

5.08 

-212,7 

88.0 

1 + .0 

3+073.188 

5.08 

-221,+ 

81.1 

l + .S 

30335.271 

+ .08 

-228,1 

84.2 

15.0 

2713+.20+ 

3.12 

-233,3 

87.+ 

15.5 

2++++. 038 

2.21 

-237.5 

100.5 

lb.O 

22188.585 

1.37 

-2+0.8 

103.7 

lb. 5 

20288.05+ 

.58 

-2 + 3.5 

10b. 8 

17.0 

18573.3+2 

-.13 

-2 + 5.7 

110.0 

17. S 

17288.818 

-.80 

-2+7.5 

113.1 

18.0 

15080.3+8 

-1. + 2 

-2+8.3 
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Problem No. 2: 


In this example, the feedline configuration consists of a large 
cavitation bubble situated in the line midway between the fuel tank 
outlet and the terminal impedance, It. As shown. .in Eigur.e D-3, the - 
entire line is acted upon by a structural velocity, Vjg(s). For analysis 
purposes, it is assumed that bubble acceleration effects are negligible, 
and that the bubble diameter is small compared to the total line length, 
Lj, + L, 3 . Proceeding as in Problem No. 1 with the assumption of zero 
tank exit pressure perturbation, the following system equation is 
obtained: - 


(D.2) 

Note that in this and the preceding example, the number of line com- 
ponents is the same. However, the present configuration results in 
a more complex functional representation of matrix 13: 

B = Eh D 3 C 1 + Cs = Bj. + B 2 (D.3) 


' P4 
_P 4 /Rj 


= D 3 Eh Dj, 


0 ' 

Lev 


-Vjj(s) 


Da D3 Cx -t- C3 


The following deductions can be made. 


( 1 ) 

NELM ~ 3 

( 2 ) 

JBNUM = 2 

(3) 

JTERM(l) = 3 


JTERM (2) = 1 

(4) 

K (1, 1) = 3 


K (1,2) = 2 


K (1,3) = 1 


K (2, 1) = 3 

(5) 

ITYPE(l) = 5 


ITYPE ( 2 ) = 2 


ITYPE (3) = 5 


The number of elements in the line. 

The number of Bi's that must be 
summed to form B. 

The number of matrices that comprise 
Bj. and B 2 , respectively. 

Two-dimensional integer array desig- 
nating the ordered subscripts of the 
matrices that comprise Bj. and I3 2 . 


Integers describing the type of compo- 
nent located in the Jth component 
position. 
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L<i-a 



L, = 15 ft. 


L 3 = 15 ft. 


R 


“|vj(s) 


o 


m 


2 

3 


} Vj( ( s ) 


Component No. Component Matrix Equations 


CD 




© 


Mean Flow Velocity 
Bubble Radius 
Line Radius 


= 50 ft. /sec. 
= 1.2 in., 

= 4.0 in. 


Qs 




-y % (s)c 


Q 


L^3J 




= J2. 


Q 4 = P 4 7R 


^3 

Q, 


- Vjj (s) C. 


3436 


Figure D-3. Line Model For Example Problem No 0 2 



(6) BSIGN (1) = +1.0 
BSIGN (2) = +1.0 


The algebraic sign preceding the 
submatrices, Bi. 


(7) SIGN=-1.0 The algebraic sign preceding the 

excitation, V_g. 

The complete input package, followed by the program output is 
shown on the following page. The transfer function, which is the 
ratio of P 4 and V^, may be normalized with respect to the pro- 
duct of the acoustic impedance, Z c , and the line cross-sectional 
area. A, as shown in the response plot of Figure D-4. Note that, 
in this example, both line segments have the same cross-sectional 
area; however, this is not a prerequisite. 
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INPUT 


EXAMPLE PROBLEM NO. 5 

3 5 3 1 ' 5 

5 2 5 

J J 
3 2 1 

3 


+ . .f npnon+ul 

+.1 OUHOO+nj 




+ . lijonpr + ui 

+ . soooon+nn + . 

bonnoo+o? +. 

5onoon+r,n +.bnonoo+o 2 

- . i n 0 p u t > + U l 

+ , 4 b 1 0 u 0 + 0 5 +. 

220000+01 +. 

4b000n-r>2 +.347000 + 02 -.238000 + 03 

+ . 3 O 0 n 0 1 > + 0 8 

+ .bbOH(iO-Ul + . 

oonnoo+on + . 

j+nnnn+ni +.133000+08 +.000000+00 

+ .4 n o n 1 1 - n s 

+.000000+00 +. 

000000+00 +. 

224000+00 +.500000+02 

+ . 1 8 0 n Q ri + 5 

+.400000+ U 1 




+ . 12C f'Oii + u i 

+ . 4 0 0 0 C f 1 + P 1 




+ . l s i' ■ i'j r i + 1 1 2 

+ . + n o u o 0 + ' i x 




OUTPUT 








NON-DIMENSIONAL 

PHASE 

OMEGA-RAD/SEC FREQ-HZ 

TRANS-P/K 

TRANS-DB 

ANGLE-DEG 

b • 3 

1.0 

41b. b82 

-24.02 

-272.8 

3.4 

1.5 

b 2 b . 3 0 0 

-20.47 

-273.2 

12. b 

2.0 

833.82b 

-17.33 

-273.7 

15.7 

2.5 

105b. 454 

-15.34 

-274.4 

18.8 

3.0 

1277.85b 

-14.28 

-275.1 

22.0 

3.5 

1505.200 

-12.8fa 

-275.8 

25.1 

4.0 

1733.78? 

-11. bO 

-27b. b 

28.3 

4.5 

1383.077 

-10.47 

-277.5 

31.* 

5.0 

223b. 741 

-3.42 

-278.3 

34 .b 

5.5 

2502.710 

-8.44 

-273.3 

37.7 

b.O 

2783.243 

-7.52 

-280.3 

40.8 

b . 5 

3081.011 

-b . b4 

-281.4 

44.0 

7.0 

3333.212 

-5.73 

-282.5 

ft 

* 

7.5 

3741.710 

-4.35 

-283.8 

50.3 

8.0 

4113.225 

-4.13 

-285,1 

53.4 

8.5 

4513.572 

-3.31 

-28b. b 

5b. 5 

3.0 

43b?. 38? 

-2.4.3 

-288,3 

58. 7 

3.5 

54b?. 531 

-l.bb 

-230.2 

b 2 . 8 

10.0 

b023.b!3 

-.81 

-232.4 

bb.O 

10.5 

bbb8.58b 

.07 

-234.3 

bq.i 

11.0 

7402.340 

.3? 

-237.8 

72.3 

11.5 

8252.522 

1.32 

-301.3 

75.4 

12.0 

3243.42b 

2.30 

-305. b 

78.5 

12.5 

10337. Ob? 

3.33 

-310.3 

81.7 

13.0 

11718. b34 

4.3b 

-317. b 

84.8 

13.5 

131bl.338 

5.37 

-32b. 1 

88.0 

14.0 

145b4.?33 

b • 85 

-33b. 7 

31.1 

14.5 

15588.030 

7.44 

-343.5 

34.2 

15.0 

15803.520 

7. 5b 

-3.3 

37.4 

15.5 

15024.751 

7.12 

-18.3 

100.5 

lb.O 

13518.373 

b * 2 1 

-31.2 

103.7 

lb. 5 

11747. bOO 

4.33 

-41.3 

10 b. 8 

17.0 

10038.740 

3 . b2 

-50.4 

110. 0 

17.5 

8528.551 

2.20 

-57.1 

113.1 

18.0 

7233.188 

.78 

-b2 • 5 
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Terminal Impedance Amplitude 



Dimensionless Frequency , cuL/C 0 


Figure D-4. Frequency Response Of The Feedline In Example 

Problem No. 2 
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Problem No. 3: 


Setup procedures for the feedline example shown in Figure 
D-5 parallel quite closely the procedures used in Problem No. 1. 
The basic difference is in the matrix representation of the center 
component, which reflects a structural acceleration input to the 
line through a viscoelastic mounting element. The functional form 
of the system equation. 


" p 4 ■ 


* or 


= DaDgDi 


.p 4 /R_ 


-Qi- 


- a i { s ) D 3 Gg 


(D.4) 


is virtually identical to Problem No. 1. Despite this apparent re- 
dundancy, this problem affords the opportunity to incorporate a pre- 
viously unused subroutine. Subroutine EIGHT, into the analysis. 


Since most vehicles employ ^lightly damped structures, a 
rather large amplification factor, Q, was assumed in order to size 
the spring constant and damping coefficient. By definition. 



The mass term is taken to be the combined fluid mass contained in 
the horizontal limbs of the feedline. For the line dimensions shown 
in Figure D-5 and Q = 100, k and b were taken to be 88 X 10 4 
lb/ft and 45.2 lb- sec/ft, respectively. 

For this problem. 


(1) ITYPE (1) = 1; 

ITYPE (2) = 8 ; 

ITYPE (3) = 1; 

(2) SIGN = -1.0 

Input values for items (1), (2) and (5) in Problem No. 1 apply also to 
this example. The response magnitude, normalized with respect to 
LOX density and a characteristic length, is shown in Figure D- 6 . 
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INPUT 


EXAMPLE PROBLEM NO. 3 

3 I 0 1 8 

1 8 i 

? 

3 2 

+ .llJDnOU+ijl 

+.iuonoo+nj +.5onnnn+no +.buoooo+u2 +,5oonoo+oo +.booono+n2 

mnnoo+i'il + . + b I nn n + n s + .2?nooci + ni + .nunono + nr) +.34?noo+n2 - .2380110 + 03 
+ . 3 nnD"n+ity +. hboooo-ni +. 000000+00 +. 000000+00 +. 13 ^ 000+08 +. 000000+00 
+ .4n8o f JO-o5 +. nnij.ton + oo +.00 0000 + no +.000° 00+00 +.500000+02 
+ . l^onoo+ua +.4-ooooo+rji 

+ . U'OnOn+fi? +. +0 000 11 + 01 +. + 52 000 + 02 +.880000+08 
+.150000+02 +.400000+01 


OUTPUT 


NON-DIMENSIONAL PHASE 


OMEGA-RAD/SEC 

FREQ-HZ 

TRANS-P/K 

TRANS-DB 

angle-deg 

b . 3 

1.0 

22.14b 

23 ,.bO 

-3.3 

3.4 

1.5 

22 • 2b 3 

23.55 

-3.3 

12. b 

2.0 

22.453 

23.72 

-4.7 

15.7 

2.5 

22.533 

23.81 

-S.b 

18.8 

3.0 

23.007 

23.33 

-b . b 

22.0 

3.5 

23.382 

30.07 

-7.7 

25.1 

4.0 

23.82? 

30.24 

-8.7 

28.3 

4 . S 

24.347 

30.45 

-3.3 

31,4 

5.0 

24.351 

30. b4 

-11.1 

34. b 

5.5 

2 5 . b 4 b 

30.8? 

-12.4 

3?. 7 

b.O 

25.442 

31.14 

-13.8 

4-0 . 8 

b.S 

27.350 

31.43 

-15.3 

44.0 

7.0 

28.384 

31.7b 

-lb. 3 

47.1 

7.5 

23 .Sbl 

32.11 

-18. b 

50.3 

8.0 

30.300 

32.43 

-20.5 

53.4 

8.5 

32.423 

32.31 

-22. b 

5b. 5 

3.0 

34.155 

33.3b 

-24 * 3 

53.7 

3.5 

35.124 

33.85 

-27.5 

b2.S 

in.o 

38 . 3bl 

34.37 

-30.4 

b b . 0 

10.5 

40.835 

34.33 

-33.7 

bq.i 

11.0 

43.748 

35.51 

-37.5 

72.3 

11.5 

45.323 

3b. 12 

-41.8 

75.4 

12.0 

50.410 

3b. 74 

-4b. 8 

78.5 

12.5 

54.10? 

37.3b 

-52. b 

81.7 

13.0 

57.833 

37.34 

-53.3 

84.8 

13.5 

bl.30b 

38.44 

-bb.3 

88.0 

14.0 

b4 . 0 8 7 

38.83 

-75.3 

31.1 

14.5 

55.735 

33.05 

-84.3 

34.2 

15.0 

55.328 

33.08 

-33.7 

3?.4 

15.5 

54 . b2 b 

38.30 

-102.8 

100.5 

lb.O 

b2 . 03b 

38.5b 

-111.4 

103.7 

lb. 6 

58.775 

38.08 

-113.2 

10b. 8 

17.0 

55.105 

37.52 

-12b. 2 

lin.o 

17.5 

51.423 

3b. 32 

-132.2 

113.1 

18.0 

47.341 

3b. 31 

-137.4 
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R 


The line length from Station 1 to the tank exit is assumed to be 3448 

negligible. Tank exit conditions exist at Station 1. 


Figure D-5. Line Model For Example Problem No. 3 





Frequency, Hz 

Figure D-6. Frequency Response Of The Feedline In Example 

Problem No. 3 





Problem No. 4: 


In this fourth example, shown in Figure D-7, an externally- 
excited segment of line is coupled at its extremities to the remainder 
of the line through two flexible bellows. This problem was dis- 
cussed briefly in Section VI. 

Initially, the user would set up the physical problem exactly 
as shown in Figure D-7. Note that the functional form of the bellows 
equations differs only in the algebraic sign of the excitation term, 
which reflects the fact that one bellows is compressed while the 
other is expanded. The response of the turbopump inlet pressure, 

Pe , in the presence of the velocity excitation, V$, is to be deter- 
mined. The component equations are combined manually to produce 
the system equation: 


•• p 6 • 


*Zi Qq," 


= D 5 D 4 D 3 D 2 D 1 


_Pe/Zt_ 


. Qi . 


Vjj DfcC* - E^D 4 C 3 - D 5 D 4 D 3 C 2 


(D.5) 


Observe that this system equation is more complex than the system 
equations in the previous examples. The following information is 
immediately available from the system equation: 

(1) NELM = 5 

(2) JBNUM = 3; therefore JTERM (1) = 2 

JTERM (2) = 3 
JTERM (3) = 4 

and 

K (1,1) = 5 
K (1,2) = 4 
K (2,1) = 5 
K ( 2 , 2 ) = 4 
K (2,3) = 3 
K (3,1) = 5 
K (3,2) = 4 
K (3,3) = 3 
K (3,4) = 2 
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Component Matrix Equations 




, P| = Q| 



q 2 2 -u w C 2 







Stationary Line 
L = 15 ft. 
r = 4.0 in. 





, Q 6 - P 6 / ^ 


g 

Combined Turbopump - Injector - 
Engine Impedance , ^ 


10 per cent Pressure Drop Through Bellows 

3445 ' 


Figure D-7. Line Model for Example Problem No* - 



(3) BSIGN(l) = +1.0 
BSXGN (2) = -1.0 
BSIGN (3) = -1.0 

(4) ITYPE (1) = 1 
ITYPE (2) = 7 
ITYPE (3) = 5 
ITYPE (4) = 7 
ITYPE (5) = 1 

(5) SIGN = +1.0 

For computational purposes, it is assumed that Zj, = 0 and 
Zt = R. In addition, there are no dissolved gases in the propellant 
(NGAS = 0). Wall elasticity is to be reflected in the speed of sound 
calculations (NELAST = 1). The entire input package is listed on 
the following page. At this point, the reader should make the 
correspondence between each numerical input quantity and the 
FORTRAN statement governing that input. The response magni- 
tude, shown in the output listing, was normalized manually, and 
the resultant response is illustrated in Figure D-8. 
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INPUT 


example problem NO. 4 




S U 1 ! 1 

5 




, l 7 -? 

-J- 




£ ^ ‘r 

* 

S *4 3 

5 4 3 £ 





t . i ts n fi.ir + n 1 

UMinnn + ni 

loonon+oi 



+ » 1 n o n rj n + ri j 

+. soopon+n o 

+ . bonnon+o2 

+,saonno+ao 

+. bonnuo+ng 

+ . lunnon + uJt 

+ • 4 h luon + nF 

+ .?'(? D n n n + u I 

+ . uonnuo+MQ 

+ . 3 4 ? n 0 G + o g - J 

+ . 3nannU + 08 

+ . bbOMjn-ni 

+. 000000+00 

+.n000Gn+nn 

•> . 1 qpnijC + 08 + .-if i'Qij ij c U 

+ . 4nHnnp«iit; 

+ * nnijunir + an 

+ .D(JOpnc> + no 

+ . GOnnan + nn 

+ . 5 n o n n f i + n g 

+ , 1 5 0 n U n + U 2 

+ . nnnono+oi 




+ « sij'Usriij + nn 

4 . *f ii n ii 1 1 * » - ■! 5 

+.i7snno+pi 



+ • ISiHi nn + ue* 

+ . 4firj^Ofi + ni 




+ # ^f'Oniir + f^n 

+ • 4 M tj U fj 1 1 - 1 F s 

+ .1 7So>m+m 



+ • 1 t >Oni:ri + Uc' ) 

+ # if n onnsur i 





OUTPUT 


NON-DIMENSIONAL PHASE 


OMEGA-RAD/SEC 

FREQ-HZ 

TRANS-P/K 

TRANS-DB 

ANGLE-DEG 

fe.3 

1.0 

775.523 

-18. b2 

-93.5 

9.4 

1.5 

1214.54? 

-14.72 

-94.2 

12. b 

2.0 

1729.378 

-11. bS 

-95.3 

15.7 

2.5 

237b. 194 

-8.90 

-9b. 9 

18.8 

3.0 

32b4 . 015 

-b . 14 

-99.0 

22.0 

3.5 

4b44.3Sb 

-3.0? 

-102.3 

25.1 

4.0 

7257.88b 

.80 

-108.7 

28.3 

4.5 

14214.91b 

b.b4 

-127.1 

31. 4 

5.0 

222bb.0b4 

10.54 

-204.9 

•3 4.b 

5. 5 

91?b.b50 

2.84 

-249.0 

3?.? 

b.O 

4807. Sb2 

-2.7? 

-259. b 

40.8 

b.S 

2801.231 

-7.4? 

-2b4 . 1 

44.0 

7.0 

1554.485 

-12.58 

-2bb.b 

47.1 

7.5 

bll.S2? 

-20. b8 

-2b?. 2 

SO. 3 

8.0 

215. 4b2 

-29.75 

-97.4 

S3. 4 

8.5 

1025.212 

-lb. 20 

-93.7 

5b. S 

9.0 

1903. bSb 

-10.82 

-94.3 

58.7 

8.5 

2942. 7b7 

-7.04 

-95.3 

b2 • 8 

10.0 

42Rl,b99 

-3.78 

-9b. 4 

bb.O 

10.5 

bl77.8S5 

-.bO 

-97.8 

b9.1 

11.0 

9209.054 

2.8? 

-99.7 

72.3 

11.5 

15048.200 

7.14 

-102.8 

7S.4 

12.0 

31348.314 

13.51 

-110.5 

78. S 

12.5 

130410.719 

25.89 

-179.5 

81.7 

13.0 

38101.154 

15.21 

-2bl . 0 

84.8 

13.5 

50414.111 

9.79 

-2b9.4 

88.0 

14.0 

14333.930 

b . 71 

-272. b 

81.1 

14.5 

112b?. 839 

4 . b2 

-274.5 

84.2 

1S.D 

9418.832 

3.07 

-275. 9 

87.4 

15.5 

8180. 4b7 

1.84 

-27b. 9 

100.5 

lb.O 

7292.480 

.84 

-277.9 

103.7 

lb. 5 

bb24.S34 

.01 

-278,7 

10b. 8 

17.0 

b!04.130 

-.70 

-279.5 

110.0 

17.5 

5b87,711 

-1.31 

-280.3 

113.1 

18.0 

5347,520 

-1.85 

-281.0 
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Problem No. 5: 


As depicted in Figure D-9, the vertical segment of the pro- 
pellant feedline is undergoing forced changes in line length. The 
structural forces applied at Stations 2 and 3 produce-measurable 
velocities - which, in general’, differ in magnitude and/or phase. 
Normally, these velocities can be related through a transfer func- 
tion, G{s). This function has been assigned a simple complex 
variable form in the computer code: G = Gj. + iG 3 . The user may 

wish to alter this form to include a specific type of frequency 
dependence based on experimental observations of velocity excited 
lines. Such a program modification can easily be accomplished. 

The length changes imposed on the vertical line segment do 
not influence the pressure-flow relationships of the horizontal seg- 
ment. 


As in the preceding examples, the first step is to describe 
each component by the appropriate matrix expression. The matrix 
equations can then be combined to yield the overall line equation 
which, for this example, is 


" Fs ' 


-pr 


= DsDi 


JPs/R- 


-Qi- 


(D.6) 


The effect of the structural velocity at Station 3 is contained in the 
column matrix, Cp.. In computing the response of P 3 to changes 
in line length, the perturbation pressure at the exit to the fuel tank. 
Pi, has been assumed to be negligible. 

At this point, the reader should have no difficulty in con- 
structing a data input package for this problem. For example, the 
input listing may take the form shown below in which the function, G, 
was taken to be 

. • * 

f 1 ;<■ • 

v 2 \fz 

The computer input and corresponding output for this package, is 
shown on the following page. The modulus of the transfer function 
relating P 3 and V 2 is shown in Figure D-10 
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INPUT 


EXAMPLE PROBLEM NO. S 

d I 0 1 


1 s 
1 

2 

+ . Lonnoo+ni 
+. loonoo+oi 

«. lOOnQP+Ol 
+.300000+08 
+.+08000-05 
+. 100000+02 
+ . 3 ri o n o n + (i 2 


3 


+.500000+00 
+.+blOOO+U5 
+ .bb0000-0'-l 
+.000000+00 
+.+ 00000+01 
+.+ 00000+01 


+. bOOPOO+02 
+.220000+01 
+. 000000+00 
+.OPnnon+GO 

+.707000+00 


+.500000+00 

+ .ooonor.!+o n 
+.000000+00 
+ .orjnnoo + nn 

+.707000+00 


+.bOOO 00+02 
+.3+7000+02 
+.135000+08 
+.500000+02 

+.280000+00 


-.258000+03 
+ ..000000 + 00 


OUTPUT 


NON-DIMENSIONAL PHASE 


OMEGA-RAD/SEC 

FREQ-HZ 

TRANS-P/K 

TRANS-DB 

ANGLE-DEG 

b . 3 

1.0 

17.3b? 

-51.52 

-157.2 

q.+ 

1.5 

17.531 

-51. +3 

-15?.+ 

iP.b 

2.0 

17.888 

-SI. 3b 

-157.3 

15.7 

2.5 

18.150 

-SI. 2 + 

-158. S 

18.8 

3.0 

18. +2 + 

-51.10 

-153.1 

22.0 

3.5 

18.71b 

-50.37 

-155.3 

25.1 

+ .0 

13.030 

-50.82 

-IbO.B 

28.3 

+ .5 

13 . 3b3 

-50. b7 

-151.7 

31. + 

5.0 

13.738 

-50.51 

-182.8 

3 + .b 

5.5 

20.1 + 0 
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Figure D-10. Frequency Response Of The Feedline In Example 

Problem No. 5 
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